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This  dissertation  describes  a physically  based  semi- 
numerical  static  and  dynamic  model  for  nonplanar-drift 
lateral  double-diffused  MOS  transistor  and  the  development  of 
a predictive  circuit  simulator  for  power  IC  TCAD.  With  the 
support  of  PISCES  simulations  for  device  operation,  the  model 
is  based  on  a regional  approach  and  its  implicit  equations  are 
solved  numerically.  A simplified  quasi -two -dimensional  analy- 
sis is  used  to  characterize  the  nonplanar  drift  region  which 
is  a feature  in  the  investigated  LDMOST  structure.  Transient 
behaviors  are  modeled  with  quasi-static  charge-based  circuit 
elements,  including  the  graded-channel  charge  dynamics  for 
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the  first  time.  The  complete  model  has  an  ability  to  account 
for  the  unique  LDMOST  characteristics  such  as  the  doping- 
graded  channel,  the  nonplanar-drif t region,  the  space-charge- 
limited  current  flow  in  the  drift  region,  and  the  inherent 
parasitic  BJT  action.  The  validity  of  the  models  has  been 
evaluated  in  the  calculation  of  device  currents,  conduc- 
tances, charges,  and  capacitances.  With  a set  of  40  model 
parameters  evaluated  in  a systematic  way,  the  computed  results 
are  in  a satisfactory  agreement  with  experimental  measure- 
ments. The  developed  model  has  been  iir^jlemented  directly  in 
SPICE2G.6  source  code  by  adding  3 new  subroutines  for  LDMOST 
model  and  modifying  18  existing  subroutines  to  accept  the  new 
model.  Unique  implementation  problems,  e.g.,  techniques 
necessary  to  ensure  current /charge  continuity,  are  described. 
The  LDMOST  is  listed  with  ID=15  in  SPICE  list  contents  with 
Y-prefix  element  name  in  the  device  card,  and  the  model  param- 
eters (ID=25)  are  listed  with  the  device  type  NLDMT  in  the 
model  card.  The  utility  of  the  new  model  in  SPICE  is  demon- 
strated through  representative  circuit  simulations,  and  the 
theoretical  predictions  are  in  qualitative  agreement  with 
other  experimental  investigations.  With  an  advantage  of 
directly  using  physical/technological  parameters  in  its 
device  model,  the  modified  SPICE  can  be  effectively  used  in 
computer-aided  power  IC  design.  The  developed  LDMOST  model  is 
also  portable  to  other  SPICE-like  circuit  simulators. 


CHAPTER  1 
INTRODUCTION 


Due  to  various  technological  advances  in  VLSI  and  high- 
voltage  devices,  power  integrated  circuits  (PIC's)  or  power 
ASICs,  which  combine  integratable  high-power  output  devices 
with  low-voltage  logic  and  control  elements  on  the  same  chip, 
have  received  increased  attention  in  recent  years.  A commonly 
used  high-power  component  in  these  circuits  is  the  lateral 
double-diffused  MOS  transistor  (LDMOST)  because  of  its  ease 
of  integration  with  CMOS  circuitry.  Furthermore,  the  perfor- 
mance characteristics  of  LDMOST' s are  generally  superior  to 
those  of  power  bipolar  transistors:  significantly  faster 
switching  time,  simpler  drive  circuitry,  the  absence  of  a 
second-breakdown  failure  mechanism,  the  ability  to  be  paral- 
lelled, and  stable  gain  and  response  time  over  a wide 
temperature  range.  Therefore,  the  availability  and  maturity 
of  LDMOST  have  made  the  PIC's  viable  and  attractive  in  many 
applications  such  as  high-voltage  displays,  motor  control- 
lers, regulators,  TV  circuits,  telecommunication,  automobile, 
switches,  and  interface  chips. 

The  technology  computer-aided  design  (TCAD)  of  these 
PIC's  requires  the  circuit-level  LDMOST  model  which  is  imple- 
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merited  into  a circuit  simulator  like  SPICE.  However,  despite 
the  impressive  progress  made  in  the  power  device  and  IC 
technology,  computer  simulation  tools,  which  are  essential  to 
the  design  and  optimization  of  power  circuits,  have  not  kept 
up  the  pace.  Although  power  devices  are  currently  modeled 
using  combinations  of  conventional  SPICE  models  with  passive 
elements,  they  are  inadequate,  and  even  with  their  empiricism 
cannot  adequately  account  for  the  unique  behaviors  of  power 
devices.  Furthermore,  the  main  disadvantage  of  these  sub- 
circuit models  is  that  they  require  fitting  parameters  which 
cannot  relate  to  device/process  parameters.  What  is  needed  is 
a physically  based  predictive  SPICE-conpatible  model,  so  that 
the  critical  model  parameters  depend  on  layout,  dopant 
profiles  and  other  device  information.  Transient  behavior 
should  be  modeled  using  charge-based  circuit  elements  rather 
than  capacitances.  This  ensures  charge  conservation  in 
circuit  simulation.  The  expressions  for  the  current  and  charge 
should  be  continuous  for  all  practical  bias  conditions  to 
avoid  convergence  problems  in  circuit  simulation.  Such  models 
can  reliably  predict  the  static  and  dynamic  performance  of 
PIC's  over  a wide  range  of  operating  conditions. 

Several  circuit-level  LDMOST  models  have  been  reported 
for  various  structures  in  the  literatures.  The  models  [Kui92, 
Poc76]  assume  a conventional  low-voltage  MOS  model  for  graded 
channel  region,  and  the  voltage  drop  in  drift  region  is 


3 


accounted  for  by  a siirple  empirical  parameter  or  a fixed 
resistor.  These  models  are  inaccurate  for  contenporary  power 
IC  simulation.  Claessen  and  van  der  Zee  [Cla86]  reported  a 
LDMOST  model,  taking  into  account  the  effect  of  non-homoge- 
neous  doping  in  the  channel  and  modulation  of  the  drift  region 
resistance.  But,  their  model  is  limited  to  dc  conditions  and 
it  does  not  account  for  the  inherent  BJT  action.  Another 
model,  which  describes  the  lightly-doped  drain  LDMOST  struc- 
ture, has  been  reported  by  Kim  et  al.  [Kim91] , in  which  the 
equation  for  the  channel  current  saturation  was  derived  based 
on  the  electron  velocity  saturation  at  the  source  of  graded 
channel.  But,  this  equation  gives  some  deficient  results  ( 
e.g.,  a constant  trans conductance  for  saturation  currents  ). 
They  used  the  conventional  bipolar  model  ( integral  charge- 
control  relation  ) to  describe  parasitic  BJT  effects,  which 
is  invalid  for  low  current  gain  (P)  devices  like  the  parasitic 
BJT  in  the  LDMOST' s.  Also  they  neglected  the  charges  associ- 
ated with  intrinsic  channel  in  their  dynamic  model. 

In  this  study,  a sophisticated  model  for  a new  LDMOST 
structure,  i.e.,  nonplanar-drif t lateral  double-diffused  MOS 
transistor  (ND-LDMOST)  shown  in  Fig.  1.1,  is  presented.  The 
purpose  of  this  dissertation  is  not  only  to  develop  a predic- 
tive device  model  accounting  for  the  unique  features  of  ND- 
LDMOST,  but  also  to  inplement  the  developed  model  into  a 
SPICE-like  circuit  simulator.  This  work  will  provide  the 


4 


SUB 


SUB 


Fig.  1.1  Cross  section  of  the  investigated  n-channel  LDMOST 

structure  showing  the  four  terminals:  the  drain 
(D) , the  gate  (G) , the  source  (S)  , and  the  substrate 
(SUB) . The  definitions  of  terminal  current  direc- 
tion and  device  geometric  dimension  are  indicated. 
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prediction  of  device  performance  as  well  as  the  computer-aided 
analysis  of  LDMOST  integrated  circiuts.  The  major  contribu- 
tions made  in  this  work  are  as  follows: 

(1)  the  modeling  of  a new  LDMOST  structure,  especially  the 
nonplanar-drift  region  by  quasi-two-dimensional 
analysis ; 

(2)  the  refinement  of  the  previous  work  for  the  graded- 
channel  model  to  properly  characterize  the  channel 
current  behavior; 

(3)  the  numerical  linkage  of  the  two  regional  ( the  channel 
and  the  drift  ) models; 

(4)  the  modeling  of  the  intrinsic  channel  charges  accounting 
for  the  graded-channel  charge  dynamics; 

(5)  the  experimental  support  of  the  developed  models  with 
measurement  data;  and 

(6)  the  implementation  of  the  developed  ND-LDMOST  model  in 
SPICE2G.6  source  code,  including  demonstration  of  the 
utility  of  the  modified  SPICE  through  representative 
power  circuit  simulations. 

In  Chapter  2,  a semi-numerical  static  model  is  described. 
The  device  is  divided  into  two  regions  ( the  graded-channel 
and  the  nonplanar-drift  ) and  the  carrier  transport  equations 
for  each  region  are  formulated.  Then,  they  are  linked  together 
by  a numerical  technique  to  provide  the  complete  solution  for 
the  current  at  a given  bias . Although  the  followed  approach 
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for  modeling  the  graded- channel  is  based  on  the  previous  work, 
the  model  has  been  refined  to  properly  characterize  the 
channel  current  behavior.  The  model  equations  for  the  nonpla- 
nar-drift  region  are  developed  based  on  physical  insight  of 
device  operation  by  a two-dimensional  (2-D)  device  simulator 
( PISCES-2B  [TMA90]  ) . The  structure  is  quite  complex,  but  a 
simplified  quasi -two -dimensional  carrier  transport  equation 
is  derived  and  solved  by  a numerical  method.  With  model  param- 
eters estimated  in  a systematic  way,  the  theoretical 
predictions  of  the  model  are  compared  to  the  experimental 
data,  and  found  to  be  in  satisfactory  agreement  over  a wide 
range  of  device  operating  conditions. 

Chapter  3 describes  the  parasitic  vertical  pnp  BJT 
action,  inherent  to  this  LDMOST  structure.  Using  the  high- 
injection  ambipolar  transport  equation,  the  parasitic  model 
is  described  and  incorporated  into  the  whole  device  static 
model.  In  addition,  the  dyneimic  properties  are  represented  by 
a quasi-static  model  of  the  integrated  charge  stored  in  the 
device.  Accounting  for  the  structural  uniqueness,  closed-form 
expressions  for  the  regional  charges  are  formulated  in  terms 
of  terminal  voltages  and  used  to  calculate  the  capacitances . 

In  Chapter  4,  the  developed  ND-LDMOST  model  is  imple- 
mented into  the  commercially  available  SPICE2G.6  program  to 
allow  an  accurate  simulation  of  lateral  DMOS  transistor 
circuits.  The  procedures  used  to  incorporate  the  model  into 
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SPICE2G.6  and  various  techniques  necessary  to  ensure  conver- 
gence are  described.  Within  the  modified  SPICE,  the  static  and 
dynamic  behavior  of  the  ND-LDMOST  are  completely  described  by 
a set  of  40  model  parameters.  Exairples  of  representative  power 
circuit  simulations  with  the  new  model  agree  well  qualita- 
tively with  theoretical  prediction  and  other  experimental 
investigations . 

In  Chapter  5,  the  main  acconqplishments  of  this  disserta- 
tion are  summarized  and  areas  for  future  work  are  suggested. 

Appendices  A and  B detail  the  linked-list  specifications 
and  nodal  admittance  matrix  of  lateral  DMOS  transistor.  Appen- 
dix C lists  a complete  source  code  of  model  routine  for 
lateral  DMOS  transistor  in  the  modified  SPICE. 


CHAPTER  2 

SEMI -NUMERICAL  STATIC  MODEL 
2.1  Introduction 

Metal -Oxide-Semiconductor  (MOS)  devices  have  been  widely 
used  for  over  two  decades  in  low  power  applications, 
especially  in  very/ultra  large-scale  integration  (VLSI/ULSI) 
circuits . Accurate  models  for  these  low  power  MOS  devices  are 
widely  referenced  in  the  literatures  [She87,  Lee88,  Yu88, 
Boo91] . However,  these  models  are  unsuitable  for  high-voltage 
double-diffused  MOS  transistors  (DMOSTs) . Two  major  changes 
in  the  conventonal  low-voltage  MOSFET  structure  are  responsi- 
ble for  the  evolution  of  the  power  MOSFET.  One  is  the  use  of 
self-aligned  double -diffused  techniques  to  achieve  a very 
short  channel  length,  which  allows  high  channel  packing 
density,  resulting  in  high  current  capability  and  low  on- 
resistance.  The  other  is  the  incorporation  of  a lightly-doped 
region  between  the  channel  and  the  drain  allowing  high  block- 
ing voltages . 

Attempts  to  model  the  DMOST's  have  been  complicated  by 
the  existence  of  these  major  features  in  structure:  the 
dopant-graded  channel  and  the  lightly-doped  drift  region. 
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While  several  physics-based  investigations  have  been 
addressed  [Dar86,  Par92,  Dol92,  Liu93],  a comprehensive  and 
practical  method  for  modeling  power  MOS  transistors  has 
remained  elusive  yet. 

This  chapter  presents  a physical  and  semi-nximerical 
model  for  the  characterization  of  power  MOS  transistors.  The 
methodology  is  based  on  a regional  device  model  concept 
[Jeo89]  that  may  have  both  theoretical  and  practical  advan- 
tages for  physical  modeling  of  complex  device  structures.  The 
investigated  ND-LDMOST  is  divided  into  two  regions  { the 
graded-channel  and  the  nonplanar-drift  ) and  the  carrier 
transport  equations  for  each  region  are  formulated.  Then,  they 
are  linked  together  by  a nxunerical  technique  to  provide  the 
complete  solution  for  current  at  a given  bias.  In  Section  2.2, 
the  graded-channel  model  is  described.  Although  the  followed 
approach  is  based  on  the  previous  work  [Kim90],  the  model  is 
refined  to  properly  characterize  the  channel  current  behav- 
ior. In  Section  2.3,  the  model  equations  for  the  nonplanar- 
drift  region  are  developed.  The  structure  is  quite  complex. 
However,  based  on  PISCES  simulations,  a siirplified  quasi-two- 
dimensional  carrier  transport  equation  is  derived  and  solved 
by  a numerical  method.  The  achieved  new  model  explains  many 
unique  properties  of  power  MOSFET's  using  arguments  based  on 
familiar  MOS  theory  and  the  physical  structure  of  the  device. 
It  accurately  represents  the  static  performance  of  power 
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MOSFET's;  this  is  confirmed  by  good  agreement  between  simula- 
tion and  experiment  results. 

2.2  Graded-Channel  Analysis 

This  section  outlines  the  formulation  of  the  static 
channel -current  equation.  For  Vjjg>0,  three  regions  of 
transistor  operation  are  defined:  the  linear  region  ( Vgg>Vrp 
and  V^h^^ch(sat)  saturation  region  ( VQg>V,p  and 
^ch^^ch(sat)  cutoff  region  ( Vgg^V,p;  in  this  case, 
the  current  conduction  through  the  channel  is  simply  zero  ) . 

Here,  V^s  and  Vqs  ane  the  applied  drain- to-source  and  gate-to- 
source  voltage,  respectively,  and  is  the  channel  threshold 
voltage.  V^h  is  the  local  voltage  drop  between  the  channel  end 
( X = Lch  ) and  the  source  ( x = 0 ),  shown  in  Fig.  2.1,  which 
is  dependent  on  the  applied  biases.  Vch(sat)  defines  the 
transition  from  the  linear  to  the  saturation  region.  Although 
the  subthreshold  region  might  be  defined,  it  is  not  important 
for  high-voltage,  high-current  devices  like  DMOST's. 

2.2.1  Linear  Region 

The  channel  current  equation  for  the  linear  region  of 
[Kim90]  has  been  refined  and  used  in  this  study.  The  model 
was  derived  based  on  five  basic  assiomptions : 1)  the  later- 
ally nonuniform  channel  (p-body)  doping  density  can  be 
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Fig.  2.1  Schematic  representation  of 
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approximated  as  a simple  exponential  function 

Nab(x)  = N^o®^P(-^ch^/Lch)  (O^x^L^^)  where  ( 
= In  (Nao/Nj3J))  ) defines  the  doping  gradient  of  channel 
region,  and  N^o  and  N^d  are  the  peak  doping  density  at  the 
source  and  the  drift  region  doping  density,  respectively;  2) 
the  strong  inversion  in  the  entire  channel  region;  3)  the 
gradual  channel  approximation  (GCA)  is  valid;  4)  the  surface 
potential  of  graded- channel  varies  linearly  with  position; 
and  5)  the  simple  relationship  between  the  inversion  charge 
density  Qn  and  the  local  potential  V;  that  is. 


(2.1) 


where  Cqx  (=  gate  oxide  capacitance  per  unit 
area,  is  the  permittivity  of  oxide,  and  is  the  average 
body  depletion  capacitance  in  the  channel. 

Although  the  model  [Kim90]  apparently  works,  it  has  two 
shortcomings.  First,  it  provides  the  abrupt  transition  from 
the  linear  into  saturation  region  by  using  a piecewise-con- 
tinuous  electron  velocity  expression.  Second,  it  overesti- 
mates Cjj  by  equating  it  to  the  value  of  body  depletion 
capacitance  at  the  source.  Instead,  in  this  study,  we  adopt  a 
continuous  electron-velocity  expression  [Sod84] , and  is 
taken  as  a fractional  value  at  the  source: 
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Cd  = Cdo[i-^^^^^]  (2-21 

where  CpQskj3^qegN^Q/(4j<|>g(0)|)  , kf3  is  a pareuneter  for  the 
depletion  capacitance,  the  body  Fermi -potential 

|4>b(^)|  = (kT/q)  In  [N^g  (x)/nj^]  , q is  the  electronic  charge,  tg 
is  the  permittivity  of  silicon,  kT/q  is  the  thermal  voltage, 
and  n^  is  the  intrinsic  carrier  density. 

From  the  above  refinements,  the  channel  current  for 
linear  region  can  be  written  by 


^^*”ox(^GS  ■*‘^DO^ch|^B  2 ^^0^  ^chj^ch  (2.3) 


where 


nch 


^^nchO 

1 + 0(Vgs-Vt) 


Vt  - ^PBch 


OX 


(2.4) 


(2.5) 


In  Eqs . (2.3) -{2. 5),  W2  is  the  device  width,  IX^chO  low- 
field  channel  electron  mobility,  and  0 are  fitting  param- 
eters for  the  effective  electron  mobility,  Vpgch  is  the  channel 
flatband  voltage,  ) is  the  critical  electric 
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field  in  channel,  and  Vgat  is  the  electron  saturation 
velocity. 

2.2.2  Saturation  Region 

When  the  transistor  is  operating  in  saturation,  the 
lateral  conponent  of  the  electric  field  in  the  channel  is 
sufficiently  large  that  the  assumptions  of  strong  inversion 
and  GCA  are  no  longer  valid.  An  accurate  solution  of  carrier 
density  and  field  distributions  for  the  graded  channel  is 
quite  complex  and  can  only  be  obtained  from  two-dimensional 
device  simulation.  It  is  claimed  [Kim90,  Kim91]  that  the 
lateral  electric  field  is  a maximum  at  the  source,  that  is, 
the  current  saturation  is  caused  by  the  electron  velocity 
saturation  at  the  source  of  the  graded  channel.  But,  according 
to  our  simulation  in  Fig.  2.2,  the  peak  position  of  the 
lateral  field  is  dependent  on  the  applied  biases  and,  conse- 
quently, the  maximum  field  moves  along  the  channel.  This 
indicates  that  it  is  unreasonable  to  get  an  equation  for 
current  saturation  from  the  electron  velocity  saturation  at  a 
fixed  position  [Kim90,  Kim91] . 

In  this  study,  however,  saturation  is  characterized  by 
a constant  channel  current  and  a channel  voltage  at  which  the 
channel  conductance  becomes  negligible.  The  onset  voltage 
Vch(sat)  constant  channel  current  is  obtained  from 


Lateral  Electric  Field  in  Magnitude  (MV/cm) 
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Fig.  2.2  PISCES-simulated  lateral  electric  field  along  the 

graded-channel  surface. 
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CH(lin) 


dV 


ch 


0. 


^ch  ^ch(*«t) 


(2.6) 


When  the  channel  voltage  increases  above  Vch(sat)»  tran- 
sistor enters  the  saturation  region  and  the  channel  current 
is  modeled  as  Ich(sat)  evaluating  the  linear  channel  cur- 
rent at  . Applying  (2.6)  to  the  equation  (2.3), 
and  assuming  the  effective  channel  mobility  is  constant  to 
avoid  complex  calculation,  we  obtain 


^ch(sat)  ^f4 


C (V 

ox  ' ' 


GS 


Vt)  CpQT|j,jj|<|)g  (0)| 


*"ox'*'^D0 


(2.7) 


where  an  empirical  coefficient  is  introduced  (0<)Cj4^1). 
Thus,  the  saturation  channel  current  is  given  by 


X 


[Co*  (Vgs  - Vt)  + Coolchl'I'B  (“)  |1  * 


(2.8) 


which  is  a function  of  only.  This  characterization  for 
current  saturation  gives  a smooth  transition  from  the  linear 
to  the  saturation  region  and  is  well  consistent  with  experi- 
mental observation.  ( This  will  be  shown  in  Section  2.5.  ) 

Note  that  our  current  model  for  the  graded  channel 
implicitly  simulates  the  saturation  characteristics  caused 
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by  pinchoff  or  velocity  saturation  somewhere  within  the 
channel,  neither  at  source  nor  at  channel  end.  This  is  con- 
sistent with  the  PISCES  simulation  in  Fig.  2.2.  When  Tj^j^  = 0, 
Eq.  (2.6)  becomes  the  typical  pinchoff  expression,  and 
(2.3),  (2.7)  and  (2.8)  become  to  conventional  MOSFET  expres- 
sions, respectively.  Another  point  should  be  made  concerning 
the  parameter  kf3  for  depletion  capacitance.  For  continuity 
in  the  expression  for  channel  current  up  to  the  cutoff 
region,  we  approximate  kf3,  not  as  a constant,  but  as 


Ic  = Ic 
^f3  ^fO 


lOk 


( for  Vqs^I.IV^  ) 


(2.9) 


fO 


- V^) 


Y ’ GS  ’T' 
Vt 


( for  V^^VQg<l.lV^  ) 


where  a constant  kfg  is  limited  to  0<  < 1 • This  Vqs -depen - 


fO 


GS 


dent  kf3  near  the  threshold  voltage  ensures  Ich  = icH(sat)  = ^ 
at.  Vqs  = V^,  which  is  important  for  the  continuity  of  the 
channel  current  for  all  operating  conditions. 


2.3  Nonplanar-Drift  Region  Analysis 


The  voltage  sustaining  capability  of  the  LDMOS  transis- 
tor is  critically  dependent  on  the  drift  region  structure.  In 
this  study,  the  investigated  device  has  a nonplanar  drift 
structure,  and  the  extended  gate  over  the  thick  field  oxide 


controls  the  direction  of  carrier  flow.  PISCES  simulations  in 
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Figs.  2.3(a)  and  (c)  show  the  current  flow  in  the  ND-LDMOST. 

In  the  laterally  graded-channel  region,  the  electron  current 

is  confined  at  the  silicon  surface.  In  the  drift  region,  the 

current  flow  is  two-dimensional,  and  the  conduction  channel 

area  varies  with  the  applied  biases.  From  Fig.  2.3,  we  see 

that  the  drift  region  can  be  divided  into  two  sub-regions:  a 

space  charge  region  (SCR)  and  a quasi -neutral  region  (QNR) . 

When  Vpg ^ V^g - Vpg^j.  ( see  Figs.  2.3(a)  and  (b)  ),  where  Vpg^jj. 

is  the  (effective)  flatband  voltage  in  drift  region,  the 

surface  area  of  the  drift  region  under  the  poly-gate  is 

partially  or  fully  depleted,  and  the  boundary  between  the  SCR 

and  the  QNR,  as  indicated  in  Fig.  2.3(b),  is  an  equipotential 

line  ( *Vgg).  When  V^g < V^g - ( see  Figs.  2.3(c)  and  (d) 

) , although  the  electric  field  due  to  the  depleted  space- 

charge  of  the  reverse-biased  p-body/n-drift  junction  has  a 

tendency  to  push  carriers  outside  of  SCR,  the  much  greater 

vertical  field  by  the  gate  accumulates  electrons  at  the 

surface  area  below  the  oxide,  and  most  of  the  voltage  drop  in 

the  drift  region  appears  predominantly  across  the  SCR,  as 

shown  in  Fig.  2.3(d).  When  the  LDMOSTs  are  merged  with  low- 

voltage  CMOS  devices,  the  epi-layer  for  drift  region  has  a 

16  _3 

relatively  high  doping  density  ( typically  - 10  cm  ) because 
it  also  serves  as  the  bulk  in  p-MOS  transistors  [Wil86,  Par88]  . 
Therefore,  the  low  current  drivability  of  LDMOST  induces  a 
negligible  voltage  drop  in  the  QNR. 
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PISCES-simulated  two-dimensional  current  flow  lines 
( (a)  and  (c)  ) and  potential  contours  ( (b)  and  (d) 
) in  the  ND-LDMOST;  (a)  and  (b)  for  V^s  = 20  V,  V^g 
= 10  V,  Vgg  = -5  V,  and  (c)  and  (d)  for  Vpg  = 8 V, 
Vqs  = 10  V,  Vgg  = -5  V,  respectively,  where  Vgg  is 
the  substrate-to-source  voltage.  In  equipotential 
contours,  the  voltage  interval  is  1 V. 


Fig.  2.3 
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Fig.  2.3 


continued 
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To  model  this  complex  voltage  drop  in  the  nonplanar  drift 
region  ( neglecting  the  voltage  drop  in  the  QNR  ) , we  define 
the  Q-point  which  gives  an  effective  two-dimensional  current 
path  in  the  SCR.  The  Q-point  is  a physical  point  on  the  bound- 
ary line  between  the  SCR  and  QNR,  and  makes  the  one- 
dimensional length  from  the  channel  end  ( x = L^h  and  y = 0 
in  Fig.  2.1  ) to  the  Q-point  become  the  shortest  route  for 
electrons  to  penetrate  through  the  SCR.  Therefore,  the  voltage 
drop  in  SCR  appears  between  these  two  points.  Obviously,  the 
Q-point  is  moving  around  the  drift  region  depending  on  the 
gate  and  drain  applied  biases.  Once  finding  the  Q-point,  we 
obtain  the  voltage  drop  in  the  drift  region  (Vqr)  by  solving 
a system  of  one-dimensional  (1-D)  differential  equations  for 
electric  field  and  potential  simultaneously. 

2.3.1  Effective  2-D  Current  Path 

A schematic  representation  of  the  nonplanar  drift 
region  is  shown  in  Fig.  2.4.  In  its  representation,  we  assume 
that  vertical  and  horizontal  diffusion  lengths  are  same  both 
in  n'^-source  and  p-body,  and  that  the  thick  oxide  is  linearly 
tapered  to  the  thin  gate  oxide.  Close  examination  indicates 
that  one  curve  and  five  lines  have  influence  on  the  electron 
transport  for  ^ V^g  - , as  indicated  in  Fig.  2.4.  A 
curve  and  three  lines  are  depletion  boundary  edges,  and  the 
remaining  two  lines  are  virtual  lines  dependent  on  dimen- 
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Fig.  2.4  Schematic  representation  of  the  nonplanar  drift 

region  in  the  ND-LDMOST.  The  specified  depletion 
boundaries  are  for  V^g  ^ V^g  - Vpg^j. . The  definitions 
for  geometric  dimensions,  curve,  and  lines  are 
indicated. 
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sions  of  the  drift  region.  Curve-1  is  the  depletion  boundary- 
edge  into  drift  region  due  to  p-body/n-drift  junction.  Line- 
2 and  line-4  are  the  depletion  boundaries  under  thin  oxide 
and  thick  oxide,  respectively.  Line-3  is  assiomed  so  that  the 
depletion  edge  under  tapered  oxide  varies  linearly  from 

to  Yfox-  When  Yox- tox)/2 , we  assume  that  Yfox  = ® 
the  depletion  boundarY  under  thin  oxide  is  extended  into  the 
tapered  oxide.  Line-5  is  a virtual  line  connecting  two 
points:  the  channel  end  ( x = Yjb»  Y = 0 ) arid  the  lower  left 
vertex  of  thick  oxide  ( x = Xj , y = (t  fox-tox)/2  ).  Line-6 
represents  the  interface  between  tapered  oxide  and  silicon. 
In  addition  to  the  geometric  dimensions  shown  in  Fig.  2.4, 
two  more  dimensions  are  defined  for  quasi-two  dimensional 
analysis ; 


(2.10) 


(2.11) 


where  X2  is  the  coordinate  of  x-axis  for  intersection  of  line- 
4 and  5,  and  r^2  represents  a physical  1-D  length  from  origin 
to  the  lower  left  vertex  of  thick  oxide. 

From  the  1-D  depletion  approximation  ( for  curve-1, 
line-2  and  4 ) , the  six  equations  for  the  curve  and  lines  can 
finally  be  written  as 
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curve-1:  x^  + y2  = + 


line-2 : y = y 


ox 


line-3:  Y = Yox"  (x-b^)  { for  Yox>  tox>/2  ) 


= y 


OX 


( yox^  (tfox-tox)/2) 


line  4.  y yfox"^  ^^fox 
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2.20) 
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(2.24) 


where  is  the  built-in  potential  of  body-drift  junction 
and  Cfox  i®  thick  oxide  capacitance  per  unit 
area.  Notice  that,  to  prevent  discontinuity  in  the  vertical 
movement  of  line-4  according  to  the  applied  biases,  we  use  a 
modifided  yfox  expression  of  (2.20)  which  results  from 
subtracting  a small  finite  depletion  thickness  under  the  thick 
oxide  to  make  yfox  = 0 when  yox  = (tfox“ 

Yds  < Vgg  - VpBdr ' ^ox  yfox  become  zero,  and  only  one  curve 
and  three  lines  ( 1,  2,  5 and  6 ) affect  the  electron  transport 
where  line-2  is  fixed  at  the  interface.  These  curve  and  lines 
determine  the  two-dimensional  path  for  carrier  transport 
through  the  SCR. 

Depending  on  the  constraint  boundary  to  be  satisfied. 
Table  2. 1-2. 3 illustrate  the  Q-point  for  eight  possible  cases 

when  Vds -^Gs  “^FBdr  three  cases  when  V^g  < Vgg -Ypp^r  • t;he 
Tables,  the  dotted  lines  ( obtained  from  curve-1  and  line-2, 
3,  and  4 ) are  approximate  boundary  edges  between  SCR  and  QNR, 
and  the  arrows  show  the  effective  2-D  current  path  in  SCR's. 
Xq  and  Yq  are  the  coordinate  of  Q-point  in  the  x-y  plane.  In 
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TABLE  2.1 

SCHEMATIC  ILLUSTRATIONS  AND  COORDINATES  OF  Q- POINT 
( For  Vpg  ^ V^g  - Vpg^j.  and  Vox^  (^fox~  ^ 
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TABLE  2.2 

SCHEMATIC  ILLUSTRATIONS  AND  COORDINATES  OF  Q-POINT 
( For  and  Yox>  (^fox"  tox)/2  ) 
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TABLE  2.3 

SCHEMATIC  ILLUSTRATIONS  AND  COORDINATES  OF  Q- POINT 

( For  Vjjg  < Vgg -Vpg^j.  ) 
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each  case,  the  Q-point  ( the  tip  of  arrow  in  the  Tables  ) is 
obtained  from  an  intersection  of  curve-1  and  one  of  the  other 
lines,  accordingly  1-D  length  from  channel  end  to  SCR  boundary 
edge  is  minim\am.  A given  V^s  and  Vqs  employs,  by  the  constraint 
boundary,  the  only  one  Q-point  among  described  eleven  cases. 
Although  in  Table  2.2,  the  case  for  X2>Xj  and  Xq>Xj  is  not 
expressed,  PISCES  simulations  reveal  that  the  Q-point  appears 
under  the  extended  poly-gate  for  practical  bias  conditions 
(below  breakdown) . Note  that  there  is  no  discontinuity  on  the 
two-dimensional  movement  of  Q-point  as  Vog  and  Vgg  change. 

2.3.2  Evaluation  of  Drift-Reion  Voltage  Drop 

Based  on  PISCES  simulations.  Fig.  2.5  depicts  the  drift- 
region  current  (Idr)  flowing  pattern  in  the  SCR.  The  principal 
2-D  current  flows  from  the  end  of  the  channel  to  the  Q-point, 
where  it  transports  the  effective  1-D  length  of  Wg^^  { = 
7(XQ-Yjb)^  + yJ  ) with  an  angle  P { = cos”^  [ (Xg-yj^) /WgcR]  ). 
The  current  spreads  out  on  both  sides.  The  lower  side,  in 
which  nearly  a half  of  Ig^  traverses,  spreads  with  an  angle  y 
( «30°)  which  is  assxxmed  to  be  independent  of  the  bias  condi- 
tion and  the  geometrical  dimension  variation  in  accord  with 
the  PISCES  simulations.  Whereas,  the  spreading  angle  of 
current  conduction  for  the  upper  side,  in  which  the  other  half 
of  IpR  traverses,  varies  with  the  applied  biases.  When  P = 0 
( refer  to  Fig.  2.3(c)  ),  the  upper-side  half  current 
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Fig.  2.5  The  SCR  current  flowing  pattern  for  the  one  half  of 

Iqpj  (lower  side)  represented  in  the  u-v  coordinate 
system. 
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traverses  with  infinitesimal  thickness,  confined  at  the 
silicon  surface.  When  both  sides  have  almost  the  same 

form,  therefore  the  shape  of  total  current  flowing  looks  like 
an  isosceles  triangle  ( refer  to  Fig.  2.3(a)  ). 

Focusing  only  on  the  lower  side,  the  mobile  sheet 
charge  density  Qn(u)  along  the  principal  current  path  can  be 
expressed  as 

Qn(u)  = -q(ho  + utany)n(u)  . (2.25) 

In  (2.25),  the  average  electron  density  n(u)  is  defined  by 

.(ho  + utanY) 

I n (u,  v)  dv 

ii(u)5-2 (2.26) 

ho  + utany 

and  the  initial  conduction  channel  thickness  ho  is  assumed  to 
be  related  to  the  Debye  length  ( k^j  JtgkT/ (q^Noo)  ) • Neglect- 
ing the  diffusion  current,  then  the  classical  carrier  trans- 
port equation  relating  the  mobile  charge  density  and 
electric  field  E^j(u)  yields 

= -W,Q„(u)  — 

1 + 

where  and  ) are  the  low-field  electron 

mobility  and  the  critical  electric  field  in  drift  region, 
respectively.  Combining  (2.25)  and  (2.27)  with  1-D  Poisson's 


(2.27) 
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equation,  the  following  relationship  for  electric  field  is 
obtained: 


du 


E„(u) 


1 r b 1 + 

ho  + utany[E^(u)  ®J  Eg 


(2.28) 


where  aQ  — and  bQ  — Iqr/ . For  a 

given  current  ( equal  to  ) , the  voltage  drop  in  drift 
region  ( Vj^j^  = V(u  = Wg^,^)  - V(u  = 0)  ) is  obtained  by  solving 

(2.28)  and  the  equation  defining  the  electrostatic  potential 

^V(u)=-E^(u)  (2.29) 

simultaneously  by  a nxamerical  method  ( the  Runge-Kutta  method 
[Con80]  ).  Eqs . (2.28)  and  (2.29)  constitute  a system  of 

first-order  differential  equations  for  electric  field  and 
potential,  and  the  Runge-Kutta  (R-K)  method  attemps  to  solve 
the  equations  explicitly  with  initial  values  for  each 
variable. 

What  is  needed  is  the  boundary  conditions  for  potential 
and  electric  field  at  u = 0 . V(u=0)  is  simply  Vch*  With 

another  expression  of  channel  current 


(2.30) 


where  E^(x)  is  the  lateral  electric  field  in  the  channel,  the 
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value  for  electric  field  can  be  obtained  by  evaluating  I 


CH 


and  IpR  atu  = 0 ( i.e.,  x = L^h  and  y = 0 in  Fig.  2.1  ) . The 


ch 


currents  and  the  mobile  charge  densities  at  x = L^h  and  u = 


0,  respectively,  are  obviously  the  same,  and  thus  setting 
= =lnn(u=0)  yields 


ch 
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(2.31) 


where  the  assxomption  is  used  and  a coeffi- 
cient kf2  is  introduced.  The  equation  (2.31)  is  valid  for  any 
angle  P because  the  electron  velocity  is  continuous  at  the 
boundary  independent  of  p.  As  V^h  increases,  K(o)|  increases 
and  the  velocity  term  in  (2.31)  approaches  to 


2.4  Numerical  Linkage  of  Regional  Models 


The  channel  current  described  in  Section  2.2  is  formu- 
lated in  terms  of  the  external  bias  Vqs  and  the  channel  voltage 
drop  for  all  operating  modes;  that  is,  its  solution  is  an 
explicit  function  of  Vqs  and  V^h/ 


^CH  “ ^ • 


(2.32) 
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Whereas,  the  drift  region  analysis  in  Section  2.3  is  described 


in  terms  of  the  drift-region  current  ( which  is  equal  to 


icH  ) • ^GS'  ^DS'  and  Vch? 


^DR  ~ ^DS  ^Ch 


~ ^ (^CH’^GS’^DS’^ch'  • 


(2.33) 


These  two  equations  (2.32)  and  (2.33)  define  a system  of  model 
equations  for  the  complete  solution  of  drain-source  current 
Ips  as  function  of  Vqs  and  V^s.  Note  that  there  are  two  unknown 
variables  V^h  and  I^h  given  terminal  biases  Vqs  and  Vqs*  We 
need  a numerical  method  to  solve  the  system  of  model  equations 
(2.32)  and  (2.33).  Their  solution  should  yield  Iqs»  nnd  a 
regional  voltage  V^h  which  enables  calculations  of  the  intrin- 
sic channel  charges  ( to  be  discussed  in  Chapter  3 ) . 

The  basic  algorithm  is  given  in  Fig.  2.6.  For  normal 
operating  bias  conditions  ( Vpg^O  and  Vgg>V^  ),  an  initial 
guess  is  made  for  the  drift  region  voltage  drop  Vqr,  and  the 
channel  voltage  drop  is  set  as  V^^  = V^g-Vj^j^  for  a given  Vqs* 
Depending  on  the  value  of  V^,jj  with  respect  to  which  is 
determined  exclusively  by  Vqs»  the  channel  current  is 
completely  defined  as  either  linear  region  or  saturation 
region.  The  calculated  graded-channel  current  is  always  equal 
to  the  drift  region  current,  consequently  the  drain-source 
current  Iqs*  The  next  step  is  the  drift  region  analysis.  V qr 
is  calculated  by  the  R-K  method  for  one  of  11  possible  2-D 
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Fig.  2.6  Flow-chart  of  the  semi -numerical  model  routine  for 

the  drain-source  current  Vqr  and  V are  the 
guessed  value  and  corrputed  value,  respectively,  for 
drift  region  voltage  drop. 


36 


current  paths  determined  by  the  two  terminal  biases  Vqs  and 
Vqs.  If  the  calculated  V is  not  close  to  the  guessed  Vqj^ 
within  a specified  tolerance,  the  new  guess  for  is 
estimated  using  the  Secant  method  [McC88]  . Otherwise,  the  main 
model  solution  loop  is  exited.  For  the  cutoff  region 
(VQg<V.p),  or  for  Vpg<0  where  the  activated  parasitic  BJT 
current  overwhelms  the  reversely  driven  channel  current  ( to 
be  discussed  in  Chapter  3 ),  is  siirply  assumed  to  be  zero. 
No  channel/drift  region  anaysis  is  done  in  these  conditions. 

2.5  Model  Verification 

The  accuracy  of  the  physical  model  has  been  verified  by 
comparisons  with  measurements.  The  test  devices  validated  are 
n-channel  ND-LDMOSTs  for  intelligent  power  IC's  developed  at 
an  SRC  member  company.  The  model  parameters  related  with 
device  geometry  and  doping  were  obtained  from  layout  and 
process  information.  The  mobility  and  electrical  parameters 
have  been  evaluated  by  a systematic  way  to  fit  the  measure- 
ments . ( The  detailed  descriptions  will  be  given  in  Chapter 
4.  ) The  low-field  mobilities  were  set  differently,  in  the 
graded-channel  and  drift  region,  to  M^nchO  “ 4.0x10  cm  /V-sec 

3 p 

and  |i.  , = 1.2x10  cm  /V-sec,  respectively.  The  empirical 
coefficients  related  with  the  channel  mobility  were  selected 
to  0 = 0.197  V”^  and  = 1.5  . The  built-in  potential  of  body- 
drift  junction  was  assumed  to  = 0.7  V,  and  the  flatband 
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voltages  were  set  to  Vp^ch  = -2.57  V and  VpB^jj.  = -0.2  V,  respec- 
tively. ( VpBch  has  been  set  to  a relatively  high  value  to  fit 
the  measured  Id'Vqs  data.  This  is  because  our  threshold 
voltage  model  in  Eq.  2.5  is  inadequate  in  a real  2-D  graded- 
channel . This  can  be  cured  by  taking  directly  as  a model 
parameter,  instead  of  calculating  it  in  the  model  equation.  ) 
The  other  parameters  were  evaluated  to  kfo  = 0.1,  kfi  = 3.0, 
kf2  = 1-0,  and  kf4  = 0.94,  respectively. 

The  comparison  of  dc  output  characteristics  are 
depicted  in  Fig.  2.7(a).  The  computed  curves  give  good  over- 
all agreement  with  measurements.  Of  particular  significance 
is  the  ability  of  the  present  model  to  show  a smooth  transi- 
tion between  the  linear  and  saturation  region,  and  to  be  more 
accurately  than  the  piece-wise  channel  model  in  [Kim90, 
Kim91] . It  is  even  more  pronounced  in  Fig.  2.7(b)  by  simply 
comparing  the  two  channel  models.  Fig.  2.8  shows  the  output 
conductance  characteristics,  where  the  derivatives  of  Ip 
with  respect  to  are  also  continous  at  the  onset  of  satu- 
ration. 

A more  stringent  verification  of  the  model  is  provided 
in  Fig.  2.9,  comparing  Iq-Vcs  transconductance  character- 
istics. Since  the  transconductance,  Gpg,  denotes  the  gain  of 
the  power  MOSFET,  it  is  an  important  parcimeter  when  the 
device  is  operated  in  the  active,  or  constant  current, 
region.  Defined  as  the  ratio  of  the  change  in  drain  current 
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Fig.  2.7  Drain  current  characteristics  (a),  and  the  compar 
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Fig.  2.8  Drain  conductance  (Gq)  versus  drain  voltage  charac- 
teristics . 
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Fig.  2.9  Characteristics  of  Id“Vqs  and  transconductance 

Vqs  = 1.0  V and  (b)  V^g  = Vgg. 
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corresponding  to  a change  in  gate  voltage  ( 

the  transconductance  varies  with  operating  conditions . The 

low-field  transfer  characteristics,  taken  at  a drain  voltage 

of  1.0  V,  are  shown  in  Fig.  2.9(a).  Agreement  is  rather  poor 

> 

compared  with  the  Id~Vds  curve,  but  this  systematic  error  can 
probably  be  reduced  by  improving  the  parameter  extraction 
procedure.  The  discrepancy  in  the  near-threshold  voltage  is 
attributed  to  the  absence  of  a subthreshold  region  model  and 
the  imperfect  description  of  the  parameter  kf3  in  our  dc 
model.  However,  the  model  correctly  depicts  the  fact  that 
maximxam  transconductance  in  the  linear  region  occurs  at  low 
gate  voltage.  Obviously,  when  the  device  is  switched  fully 
on,  the  transistor  will  be  operating  in  its  ohmic  region 
where  the  gate  voltage  will  be  high.  In  that  region,  a change 
in  an  already  high  gate  voltage  will  do  little  to  increase 
the  drain  current;  therefore,  Gps  will  diminish  to  zero.  This 
feature  is  in  contrast  to  it's  behavior  in  the  saturation 
region  as  seen  in  Fig.  2.9(b)  . By  taking  V^g  = V^g  , the  value 
of  Gps  is  determined  from  the  saturation  region  of  the  Id'Vqs 
characteristics  where  a change  in  V^g  no  longer  significantly 
influences  on  Ip,.  Therefore,  only  VQg  would  make  a change  in 
the  Gpg  characteristics,  increasing  Gpg  consistently  as  V^g 
rises . 

The  on-resistance,  RDs(on)  ( = V^g/Ip,  ) , of  a power  MOS- 
FET  is  an  important  figure  of  merit  because  it  determines  the 
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amount  of  current  the  device  can  handle  without  excessive 
power  dissipation.  When  switching  the  power  MOSFET  from  off 
to  on,  the  drain-source  resistance  falls  from  a very  high 
value  to  RDs(on)'  which  is  a relatively  low  value.  The  compar- 
isons of  static  on-resistance  characteristics  are  shown  in 
Fig.  2.10.  As  the  drain  current  rises,  especially  above  the 
continuous  rating,  the  on-resistance  also  increases.  Agree- 
ment is  generally  good  except  in  the  ohmic  region  where  the 
model  values  are  about  18  % less  than  measurements.  As  Fig. 
2.10  indicates,  increasing  the  gate  voltage  has  a diminish- 
ing effect  on  lowering  the  on-resistance.  To  minimize  RDs(on)' 
the  gate  voltage  should  be  large  enough  for  a given  drain 
current  to  maintain  operation  in  the  ohmic  region.  But, 
somewhat  like  driving  a bipolar  transistor  deep  into  satura- 
tion, unnecessarily  high  gate  voltages  will  increase  turn- 
off time  because  of  the  excess  charge  stored  in  the  input 
capacitance . 

As  an  example  of  device  characteristics  at  high  temper- 
ature, Fig  2.11  shows  the  temperature  dependence  of  the 
static  on-resistance.  As  temperature  increases,  the  device 
characteristics  are  crudely  determined  by  the  variation  of 
the  thermal  voltage,  the  intrinsic  carrier  concentration  and 
the  carrier  mobility.  The  value  of  these  parameters  at  tem- 
perature T is  related  to  their  value  at  a reference  tempera- 
ture Tq  ( = 300  K ) by  the  following  equations  [Dol92] : 


RDS(on)XWz  (kohm-um) 
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Fig.  2.10  Characteristics  of  static  on-resistance. 
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Vto<T)  = (2.34) 

ni(T)  <2.35) 

^^nchO^'^)  * (2.36) 

^*„d(T)  = ^„d(To)(^)'“-  (2-37) 


In  Eqs . (2 .34)  - (2 .37)  , V^q  the  thermal  voltage,  k is  Boltz- 
mann's constant,  and  Eg  is  the  silicon  energy  gap.  Because  the 
mobility  term  dominates,  increasing  temperature  for  a fixed 
Vqs  tends  to  decrease  the  carrier  mobility,  consequently 
increases  the  on-resistance. 

2 . 6 Summary 

A predictive  static  model  for  the  nonplanar-drif t 
lateral  DMOS  transistor,  which  is  intended  for  SPICE  imple- 
mentation, has  been  described  and  verified.  A detailed 
description  of  the  model  formulation  is  presented  and  the 
assumptions  used  in  the  derivation  of  model  equations  are 
stated.  With  the  support  of  PISCES  simulations  for  device 
operation,  the  models  are  based  on  a regional  approach  and  its 
implicit  equations  are  solved  nvimerically . The  graded-channel 
model  has  been  refined  to  properly  characterize  the  channel 
current  behavior.  A simplified  quasi-two  dimensional  analysis 
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has  been  used  to  characterize  the  nonplanar  drift  region  which 
is  a feature  in  the  investigated  lateral  double-diffused  MOS 
transistor  structure.  Different  from  the  existing  power 
device  sub-circuit  models  [Sco91,  She91] , the  model  developed 
herein  has  an  ability  to  account  for  the  unique  LDMOST 
characteristics . 

The  model  has  been  tested  successfully  in  the  calculation 
of  device  currents,  conductances,  and  on-resistance,  and 
shown  continuous  behavior  under  all  operating  regions.  This 
property  of  continuity  will  improve  the  convergence  problems 
arising  in  circuit  simulation.  With  model  parameters  evalu- 
ated by  a systematic  procedure,  the  simulation  results  for  the 
static  characteristics  are  in  a good  agreement  with  experi- 
mental measurements.  The  feature  of  accuracy  is  attributed  to 
the  physical  representation  of  all  the  important  effects 
occuring  in  the  ND-LDMOST.  Because  the  model  is  physically 
based  and  the  number  of  fitting  parameters  is  chosen  as  small 
as  possible,  the  extracted  model  parameters  can  be  used  to 
monitor  and  predict  the  effect  of  process  changes.  Even  if  the 
penalty  for  computational  efficiency  due  to  the  regional 
approach  exists,  this  new  model  is  physical,  accurate,  and 
suitable  for  incorporation  into  circuit  simulators. 


CHAPTER  3 

PARASITIC  COMPONENT  AND  CHARGE-BASED  DYNAMIC  MODEL 

3.1  Introduction 

In  Chapter  2,  a semi -numerical  static  model  for  the  ND- 
LDMOST  has  been  introduced.  In  this  chapter,  an  analytical 
model  for  the  parasitic  components  in  the  ND-LDMOST  is  first 
derived,  which  is  self-consistent  with  the  semi -numerical 
model  and  essential  to  the  simulation  of  inductive-load 
switching.  Then,  the  dynamic  model  of  the  device  is  presented, 
which  describes  the  transient  current  behaviors  corresponding 
to  the  variation  in  terminal  charges.  Consequently,  the  device 
model  needed  for  computer-aided  circuit  analysis  of  the  inves- 
tigated DMOST  is  completed  by  adding  these  models  to  the 
drain-source  current  model  described  in  Chapter  2 . 

Inherent  in  most  power  MOSFET's  are  parasitic  p-n 
junctions  formed  by  the  body-to-source,  the  body-to-drif t , 
and/or  the  substrate-to-drift . In  many  applications,  these 
internal  diode  are  never  forward  biased  and  do  not  influence 
circuit  operation.  But,  depending  on  the  terminal  connection 
and/or  bias  conditions,  these  junctions  may  act  like  the 
bipolar  junction  transistor  (BJT)  during  transient  switching. 


48 


49 


causing  a significant  effect  in  circuit  operation.  Because  of 
its  extensive  junction  area,  the  current  ratings  of  the 
parasitic  BJT  are  the  Scune  as  the  MOSFET's  continuous  and 
pulsed  current  ratings.  In  Section  3.2,  a comprehensive  but 
compact  model  is  developed  for  this  parasitic  BJT,  which  is 
necessary  for  an  accurate  transient  simulation  in  power  MOSFET 
circuits . 

In  Section  3.3,  we  describe  the  development  of  a charge- 
based  transient  model  for  the  ND-LDMOST  that  accounts  for  the 
structural  uniqueness.  Based  on  the  quasi-static  approxima- 
tion, closed-form  expressions  for  the  terminal  charges, 
continuous  throughout  the  regions  of  operation,  are  derived 
in  terms  of  terminal  voltages.  The  intrinsic  graded-channel 
part  of  the  device  is  of  major  concern  here.  But,  a compre- 
hensive description  of  the  extrinsic  part  is  also  given,  and 
verified  by  comparing  simulations  with  experiment  results. 

3.2  Parasitic  BJT  Model 

There  are  two  parasitic  diodes,  which  form  a vertical  pnp 
bipolar  transistor,  in  the  junction-isolated  n-channel  LDMOST 
like  the  investigated  device:  the  body-to-drift  and  the 
substrate-to-drift  diode,  shown  in  Fig.  3.1.  When  the 
substrate  terminal  is  connected  to  the  most  negative  voltage 
in  a power  IC,  a negatively  biased  drain-to-source  voltage  Vqs 
( which  could  happen  momentarily  during  transient  switching 
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Fig.  3.1  Parasitic  components  created  by  the  body-to-drif t 

and  the  substrate-to-drif  t junction  in  the  ND- 
LDMOST . 


51 


of  the  LDMOST  ) causes  the  body  diode  to  turn  on,  so  that  the 
deep-p'''  body  becomes  a terrific  hole  emitter  for  the  activated 
pnp  bipolar  transistor  ( forward  active  mode  ) , as  shown  in 
the  PISCES  simulation  of  Fig.  3.2.  When  the  substrate  terminal 
is  connected  to  the  source  ( i.e.,  the  substrate-to-source 
voltage  Vgg  = 0 ) , the  negative  V^g  turns  on  both  diodes  and 
the  BJT  enters  into  the  saturation  mode.  Especially  for  this 
terminal  connection,  the  both  parasitic  diodes  play  an  impor- 
tant and  useful  role  as  a drain-source  clamp  diode  in 
inductive-load  switching.  In  multi-transistor  configurations, 
such  as  the  totem  pole  network,  each  transistor  is  protected 
from  excessive  flyback  voltages,  not  by  its  own  parasitic 
diode,  but  by  the  diode  of  the  opposite  transistor  [Mot89] . 
In  either  case,  the  source,  drain,  and  substrate  terminals  act 
as  the  emitter,  base,  and  collector  of  BJT,  respectively,  and 
the  parasitic  current  overwhelms  the  entire  device.  Although 
another  parasitic  transistor  ( source-body-drain  npn  BJT  ) 
might  be  turned  on  when  a large  positive  dV^g/dt  is  developed, 
the  deep-p"*"  body  diffusion  prevents  this  activation  and  thus 
the  npn  BJT  doesn't  turn  on  during  typical  switching  of  DMOST 
[Rav91] . 

The  equivalent  dc  network  representation  for  parasitic 
vertical  pnp  transistor  is  illustrated  in  Fig.  3.3.  The  ohmic 
resistance  of  deep-p"*"  body  diffusion,  drift  region,  and 
substrate  in  the  ND-LDMOST  are  modeled  by  the  linear  resistors 
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Fig.  3.2  PISCES-simulated  hole  current  vectors  in  the  ND- 

LDMOST  at  Vds  = -0.8  V,  Vqs  = 0.0  V and  Vgs  = -5.0 
V,  showing  the  parasitic  BJT  action. 
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Fig.  3.3  Equivalent  dc  circuit  representation  of  the  para- 
sitic vertical  pnp  bipolar  transistor.  The  internal 
nodes  ( emitter,  base,  and  collector  ) of  the  BJT 
are  indicated. 
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Rbd/  Rdr'  ^suB*  characteristics  of  the  intrinsic  pnp 
BJT  are  described  by  the  nonlinear  current  source  Igp  { base 
current  for  foward  active  mode  ) , IgR  ( base  current  for 
reverse  active  mode  ) , and  Iqt  ( base  transport  hole  current 
where  Icp  is  for  foward  active  and  Icr  is  for  reverse  active 
) . The  reverse  active  mode  of  operation  ( for  reverse-biased 
body  diode  and  foward-biased  substrate  diode  ) is  not  the  case 
for  practical  bias  conditions,  but  is  considered  only  for  the 
purpose  of  describing  the  operation  in  saturation. 

3.2.1  Base  Region  Analysis 

The  analysis  of  the  pnp  BJT  is  described  using  the  coor- 
dinate system  defined  in  Fig.  3.4,  where  yg,  ygR/  Yc>  Vsc 
represent  edges  at  SCR  of  the  emitter-base  (E-B)  and  collec- 
tor-base (C-B)  junction.  In  deriving  the  parasitic  BJT 
model,  basically  we  make  four  assuirptions : 1)  wide  emitter/ 
base/collector  region,  2)  high  injection  in  the  base,  3) 
superposition  of  base  mobile  carriers  like  the  Gummel-Poon 
model  [GumVO] , and  4)  negligible  generation/ recombination 
current  in  junction  SCR's. 

From  the  above  assumptions,  the  ambipolar  transport 
equation  (ATE)  in  the  base  for  steady-state  gives  the  solu- 
tion for  the  excess  hole  density; 


p(y)  = PBp(y) 


(3.1) 
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Fig.  3.4  Coordinate  system  used  to  model  the  vertical  pnp 

bipolar  transistor  in  the  ND-LDMOST. 
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(3.2) 


(3.3) 


where  L;^g  is  the  ambipolar  diffusion  length  in  the  base,  and 
W3fg  is  the  normalized  base  width  ( s (yr~yE)'^i‘AB  ^ * 

The  boundary  conditions  for  the  excess  hole  concentra- 
tion at  the  emtter  and  collector  side  are  given  by 


where  Veg  and  V^b  are  the  quasi-Fermi  potential  separation  at 
E-B  and  C-B  junction,  respectively.  In  (3.4)  and  (3.5),  the 
second  term  ( negative  1 ) in  the  brackets  is  incorporated 
for  a crude  description  of  the  current  transport  at  reverse- 
biased  junction. 

Since  the  quasi-Fermi  potential  separations  at  each 
junction  become  the  junction  voltages  when  reverse  biased, 
the  moving  boundary  yg  and  Yq  iri  the  base  can  be  approximated 
for  all  bias  regions; 


(3.4) 


(3.5) 


( for  ) 


(3.6) 


= ^jdp 


( for  VEB>Vbd  ) 
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(3.7) 


= ^epi  ( VcB^Vsub  ) 

where  Nj^g  is  the  substrate  doping  density  and  Vs„h  ^ 
= (kT/q)  In  (N^gN^p/n?)  ) is  the  built-in  potential  at  C-B 
junction . 


3.2.2  Base  Current 

The  base  current  for  each  mode  is  sinply  expressed  as 
the  sum  of  the  integrated  recombination  current  in  the  base 
and  the  back  injection  current.  For  the  forward  active  mode, 
the  low  injection  prevails  over  the  emitter  region.  Account- 
ing for  bandgap  narrowing  and  majority  carrier  degeneracy  in 
the  emitter  injection  current  iN(ysE)  [Fos81]  , it  is  given  by 


^BF  ~ 
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PBp(y)dy  + iN(ysE) 
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Tjjb  sinh(Wgj^)  L^^^WkT^'  J 


fqV 


^ expi  , 


[ 


EB 


)-] 


(3.8) 


where  Ag  is  the  effective  E-B  ( body-drift  ) junction  area, 
Xjjg  is  the  high  injection  carrier  lifetime  in  the  base, 
NAB(eff)  ( ® constant  ) is  the  effective  doping  density  in 
deep-p"^  body  region,  and  is  the  electron  lifetime  in  the 
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emitter. 

For  reverse  active  mode,  the  solution  of  another  1-D 
ATE  in  the  collector  where  high  injection  prevails,  yields 
the  excess  electron  density  ncR(y)  in  the  collector.  Then, 
the  base  current  is  given  by 
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qAcj^c 

— I PBR(y)dy  + 


Pf 

'’hc’/Ysc 


Hcr  (y)  dy 


= qA^n^ 


r^AC  1 ^AB  ^OSh  (Wgjj)  1 

ir  r^^cB^  ,1 

L^hc  "^hb  _ 

(3.9) 


where  A^  is  the  effective  C-B  ( substrate-drift  ) junction 
area,  and  and  are  the  high  injection  carrier  lifetime 
and  ambipolar  diffusion  length  in  the  collector, 
respectively . 


3.2.3  Base  Transport  Current 


Following  the  approach  of  [McD88] , the  set  of  equations 
for  base  transport  hole  current  can  be  written  as 
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(3.11) 
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where  A£^  is  the  effective  BJT  cross-sectional  area,  and  bg 
and  are  the  electron-to-hole  mobility  ratio  and  the  ambipo- 
lar  dif fusitivity  in  the  base,  respectively.  In  (3.11)- 
(3.14),  J.pp  and  define  the  emitter  current  density  for  the 
forward  active  mode  and  the  collector  current  density  for  the 
reverse  active  mode,  respectively.  Then,  the  total  base  trans- 
port hole  current  is  finally  given  by 
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1 +;~cosh(W„J 


■^EC'^RBO' 
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sinh  (W„  J 
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(3.15) 
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(3.16) 


(3.17) 


(3.18) 


This  IcT/  combined  with  Igg,  yields  the  substrate  current  in 
the  ND-LDMOST.  Notice  that  the  asymmetric  BJT  structure  (p’*’- 
n-p)  and  operation  ( i.e.,  low  injection  in  the  emitter  and 
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high  injection  in  the  collector  ) induces  the  base  transport 
current  for  V^b  = (3.15). 

3.3  Expressions  of  Regional  Charges 

The  dyncunic  properties  are  represented  by  a quasi-static 
model  in  terms  of  the  integrated  charges  stored  in  the  device. 
The  charges  in  the  ND-LDMOST  include  1)  charges  associated 
with  intrinsic  channel,  2)  charges  associated  with  the  gate 
overlap,  and  3)  charges  associated  with  junction.  Fig.  3.5 
shows  all  of  the  above  charges  associated  with  the  ND-LDMOST. 
Qgj  and  Qqi  are  channel  inversion  charges  associated  with  the 
source  and  drain  terminal,  and  Qbi  is  the  body  charge  in  the 
channel  region.  Qqsol  Qgdol  represent  the  charges  associ- 
ated with  gate  overlap  over  the  source  and  drift  region.  Qjg 
and  Qjc  are  the  depletion  charges  associated  with  the  body- 
drift  (E-B)  and  substrate-drift  (C-B)  p-n  junction,  and  Q^e 
and  Q]x;  are  the  forward-biased  injection  charges  associated 
with  those  junctions,  respectively.  The  charges,  in  general, 
are  dependent  on  node  voltages  and  the  capacitive  currents  are 
given  by  the  derivatives  of  the  charges  with  respect  to  time. 
However,  it  should  be  noted  that  the  charges  associated  with 
intrinsic  channel  region  ( Qsj,  Qdi»  Qbi  ) are  dependent 

on  Vch  and  Vqs-  V^h  is  a (virtual)  internal  node  voltage  and 
given  by  the  solution  of  Ipg  ( discussed  detail  in  Section  2.4 
) . 
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p-sub 


Fig.  3.5  Schematic  cross-sectional  view  of  the  ND-LDMOST 

with  the  associated  charges . 
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3.3.1  Intrinsic  Graded-Channel  Charges 

Based  on  the  assuitptions  and  results  in  the  steady- 
state  analysis  of  the  graded  channel,  the  intrinsic  charges 
Qsi,  Qdi,  and  Qgi  are  formulated.  Following  the  channel  inver- 
sion charge  partitioning  scheme  in  [OhSO]  , Qsi/  Qdi'  Qbi 
are  taken  as 


^si  “ 

Qdi  ” 

^BI  “ ^z 


(x)dx 


*ch 


ch 


Q^(x)dx 


I 

Jo 


(x)dx 


(3.19) 

(3.20) 


(3.21) 


where  Qd(x)  is  the  body  depletion  charge  density. 

Rearranging  the  relationship  (2.1)  to  obtain  a term  dQ^ 
and  a term  muliplied  by  dV,  integrating  these  terms  along  the 
channel,  from  0 to  x and  from  0 to  V(x),  respectively,  we 
obtain 


(3.22) 


Then,  applying  V (x)  = x , the  channel  inversion  charge 

density  is  expressed  as 
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-(Cox  + Cdo)^='1  <3-23) 

which  is  linearly  dependent  on  the  position  x,  shown  in  Fig. 
3.6(a).  The  actual  Qn(x)  distribution  in  the  graded  channel 
is  quite  complex  and  may  not  be  monotonously  increasing  or 
decreasing  with  regard  to  the  position  x.  But  this  first- 
order  characterization  of  Qn(x)  is  well  consistent  with  MOS 
physics.  For  low  Vug  ( i.e.,  low  Vch  )/  it  simulates  the  added 
inversion  charge  away  from  the  source  in  the  laterally 
graded  channel.  This  is  an  important  feature  in  the  graded 
channel,  different  from  the  conventional  MOS  channel.  As  Vch 
increases,  the  inversion  charge  in  the  channel  end  would 
decrease  because  of  high  lateral  electric  field  in  that 
position.  The  channel  enters  the  saturation  region  with  fur- 
ther increasing  Vch-  Note  that  the  pinch-off  voltage  of 
inversion  charge  at  the  channel  end  (Vch(pin))  is  always 
greater  than  Vch  (sat)  sny  Vqs-  Physically,  this  means  that 

our  first-order  model  implicitly  conforms  to  the  current 
saturation  characteristics  caused  by  pinch-off  or  velocity 
saturation  somewhere  within  the  channel  not  by  pinch-off  at 
X = Lch»  and  that  the  electron  inversion  charge  is  always 
negative  along  the  channel  for  the  linear  region.  Applying 
(3.23)  to  (3.19)  and  (3.20)  yields 

Qsi  = 
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° ^ '^chl  ^ '^ch2  ^ '^ch  (pin) 


ch2 


^ox  ^DO 


ch(pin) 


_ ^ox(  ^GS  ~ ^T ) •*•  ^^Do\h|***B 


^ox^DO 


-Qn(U) 


Vch  = V, 


ch3 


^(0) 


-Qn(Uch) 


Fig.  3.6  Inversion  charge  density  (a)  and  charge  partition- 
ing schemes  (b) , in  the  graded-channel  for  VQg>V,p. 
Note  that  ^gs>  so  that  Qn(x) 

distribution  marked  with  dashed  line  in  (a)  never 
occur  in  our  first-order  characterization  of  Q^ix.)  . 
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+ ^Coonchl'I’B («)| - (Cox  + Cm)  Veh  ] (3  • 24 ) 

«DI  = -K^'ch[SC„x(VGS-VT) 

^CDO^ch|*l*B  (0)  I “ ^ (Cox  ■*■  Cdq)  '^ch  ) (3.25) 

for  the  linear  region.  Fig.  3.6(b)  shows  the  partitions  of 
these  inversion  charges  in  the  graded-channel . 

From  the  depletion  approximation  in  the  channel,  (3.21) 
can  be  written  as 


p^ch 

« 

0 

VNab(^)J2|4)b(x)|  + ^x^ 

dx  . (3.26) 

To  find  a simple  analytic  charge  equation  in  (3.26),  the 
square  root  term  of  the  channel  doping  density  and  the  body 
Fermi -potential  ( a wealc  function  of  x ) are  approximated  as 
and  |<t>gj , where 

log(^N^gJ  = (logN^o+ logNgg)/2  (3.27) 

R = I|<l'B(0)|+|'t’B(I'ch)|l/2-  (3.28) 


respectively.  Using  these  approximations,  the  body  charge 
for  linear  region  is  expressed  as 


Qgj  = -W,L^v,K 


[ 


Vch  + 2 


z ch  BI 


V 


(3.29) 


ch 


where  Kg^  s 2,^q£gl^/3  . 
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For  the  saturation  region,  the  charges  associated  with 
the  channel  are  obtained  directly  from  (3.24),  (3.25),  and 
(3.29)  by  replacing  V^h  with  Vch(sat)  • 


(2  - k ^4)  (0)  I ] 


(3.30) 


1 


Qdi  = -6WzLch[(3-2kf4)C,^(V<,s-V^) 


+ (4  2kf4)Cj3^T)^jj|<l)g(0)|  ] 


(3.31) 


Qbi  ~ 


Wz^ch^Bl 


^ox  ^DO 


^f4  ^ox  ^DO^chl^B 


(3.32) 


X 


^ox^^GS  rolI“lV/2 


f4 


^ox  ^DO 


iq  -(^ra) 


which  are  all  dependent  only  on  Vqq. 

When  the  gate  voltage  drops  below  the  threshold  volt- 
age, the  channel  enters  the  cutoff  region.  The  inversion 
charge  is  ass\imed  to  be  negligible  compared  to  the  other 
charge  components: 

Qsi»Qdi“0-  (3.33) 


The  body  charge,  when  Vpg^jj<VQg^V.j, , is  approximated  for  a 
continuous  expression  at  all  bias  conditions  as 


Qbi  = 


r 

Ks- 


im  (3.32) 

v„ 


GS  ^FBch 


^FBch 


f 


(3.34) 
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Applying  I'Hopital's  rule  [Ste87]  to  (3.34)  to  find  a limit 
of  Qbj  in  saturation  as  Vqs  approaches  the  equation  (3.34) 
becomes 


If  the  gate  voltage  decreases  below  the  flatband  voltage, 
the  channel  area  is  acciimulated,  and  the  body  charge  is  given 
by 


3.3.2  Gate  Overlap  Charges 

The  gate  overlap  charge  over  a heavily  doped  source  is 
simply  modeled  by  a constant  capacitance: 


where  Al  is  defined  at  Fig.  1.1.  This  siitple  model  is  inad- 
equate for  the  gate  overlap  charge  over  a lightly  doped  drift 
region.  In  such  a case,  the  overlap  capacitance  varies  with 
the  applied  gate-to-drain  voltage  due  to  the  depletion  of 
mobile  carriers.  { Physically,  this  overlap  capacitance  is 
important  in  transient  characteristics  of  power  DMOS  tran- 
sistors. Because  the  capacitance  is  often  connected  between 
the  output  and  input  terminals  of  a high  gain  circuit,  it  may 


(3.35) 


(3.36) 


(3.37) 
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dominate  the  switching  speed  of  circuits  through  the  Miller 
effect.  ) 

Accordingly  in  this  study,  when  the  surface  area  of 

drift  region  is  accumulated,  Qqdol  is  expressed  as 

^GDOL  ~ t^ox^dl  ^fox^  ^fox^d3^ 

^ (^DS  ~ ^GS  ^FBdr^  ( for  Vjjg  < Vgg  — Vpg^j.  ) (3.38) 

where  the  tapered  oxide  capacitance  is  assumed  to  be 
WzLd2  + Cjq^)/2  . When  the  surface  area  is  depleted,  by  sim- 
ply using  equations  (2 .13) - (2 .15)  ( see  Fig.  2.4  ),  the  fol- 
lowing charge  equations  for  Qqdol  derived: 

®GDOL~  ~*5^QD^z^DD^ox|  ^dl  r "ZT  ^oxl 

L ""fox  ^ox  J 

{ for  . and  y ^ (t.,  -t  )/2  ) 

' DS  GS  FBdr  -'ox  fox  ox'  ' 

(3.39) 

“ “^^QD^z^DD^  (i*dl  ■*’^d2'^)y^ox‘*'  (i'd2'^'*’i'd3)yfoxl 
( for  Vgg  ^ Vgg  — Vpgi^j.  and  yQ^^  ^^fox~  ^ox^ ^ 

in  which  fQp,  ( charge  sharing  factor:  0<fgg^l)  is  intro- 

duced. This  factor  reflects  that  the  total  depleted  charges 
under  the  overlapped  gate  are  partially  associated  with  the 


gate . 
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3.3.3  Junction-Related  Charges 


The  following  equations  are  used  to  describe  Qjg  and 


Qjc  • 


QjE  ” Vbd^DD 


1- 


V 


EB 


Vbd 


= 0 


( for  VEB<Vbd  ) 


( for  ) 


(3.40) 


QjC  “ ^C.P^^sVsub 


N N 
^DD  ^AS 


1- 


V 


CB 


V 


sub  / 


{ for  VcB<Vsub  ) 


(3.41) 


= 0 


( for 


CB  - Tsub 


which  ensure  continuous  expression  at  all  bias  conditions. 

When  the  parasitic  BJT  turns  on,  the  excess  mobile 
charges  stored  in  the  device  affect  significantly  the  tran- 
sient characteristics.  These  charges,  using  the  charge  con- 
trol approximation,  are  given  by 


= PBE(y^^y+\E^n^ysE'> 

COSh(WB^)-lp  fqVEB^  ,1 
= sinh(W,J 


qA 


rqv, 


•^E^i^idpr  I '*’EB 
^AB(eff)  •-  ^ 


)->] 


(3.42) 


I 


Qpj,  = qA,,  I p„B  (y)  dy  + qA 


C I *'BR 

y, 


C I 

•'Vsc 


(y)dy 


70 


= qA^rii 


cosh  (Wg^)  - Ilf  rqVcB^  I 

sinh(Wg^)  jL^^^^2kT>'  J 


(3.43) 


from  the  base  current  expressions  of  (3.8)  and  (3.9), 
respectively . 


3.4  Model  Verification 


Simulations  have  been  carried  out  in  order  to  predict 
the  parasitic  BJT  action  and  to  verify  the  dynamic  model  in 
the  ND-LDMOST,  discussed  in  Sections  3.2  and  3.3,  respec- 
tively. Device  geometry  and  process  parameters  are  the  same 
as  those  in  Chapter  2.  Typical  diffusion  parameters  were 
assumed  that  the  electron-to-hole  mobility  ratio:  bg  = 3.5, 
the  carrier  lifetimes:  = 0.15  usee,  X^g  = 0.3  usee,  Xg,^  = 
0.3  usee,  and  the  ambipolar  diffusion  lengths:  = 8.6  um, 
i'AC  = 23.6  um,  respectively.  The  ohmic  resistances  were  esti- 
mated from  process  information:  RgQ  = 3.5  Q,  Rqr  = 56.7  , 
and  Rsub  = 62.1  , respectively. 

Fig.  3.7  shows  the  simulated  current -voltage  character- 
istics by  the  parasitic  BJT  action.  For  a given  V^s  and  Vss' 
the  respective  dc  nonlinear  nodal  equations  at  E,  B,  and  C 
nodes  in  Fig.  3.3  are  rearranged  to  find  the  internal  node 
voltages  Veb  and  Vcb^ 


gi(v 


EB’ 


(^BD  i 


DR'  BF 


4-  R T 
^ ^DR“^BR 


^BD^CT 


V + V 
''eb  ^ ''ds 


(3.44) 


(ujn/vn) 


71 


(a) 


Fig.  3.7  Simulated  parasitic  BJT  characteristics  for  the  ND- 

LDMOST:  (a)  the  substrate  current,  and  (b)  the  drain 
current . 


q/Wz  (uA/um) 
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Vds  (V) 


Fig.  3.7 


continued 


73 


92  (Veb»  “ ^DR^BF  ^SUB^  ^BR  ^SUB^CT '*’^CB  ^SS 

= 0,  (3,45) 

and  solved  by  the  Secant  method  [McC88] . ( When  implemented 

into  SPICE,  the  internal  node  voltages  are  found  in  SPICE 
nodal  analysis.  ) For  Vgg  = 0,  the  negatively  biased  Vqs 
turns  on  both  the  body-drift  and  substrate-drift  diodes, 
where  the  base  transport  hole  current  Iqt  is  negligible. 
Hence  Isub  hsve  a typical  diode  characteristics.  For 

low  negative  Vgg,  the  pnp  BJT  operates  in  saturation  mode, 
yielding  a low  substrate  current.  When  Vgg  is  sufficiently 
negative,  the  BJT  operates  in  the  foward  active  mode,  and 
IguB  doesn't  much  change  with  further  negative  Vgg.  All  of  the 
above  parasitic  BJT  behaviors  are  well  depicted  in  Fig.  3.7. 

For  the  validation  of  the  dynamic  model,  normalized  plots 
of  the  four  intrinsic  channel  charges  are  shown  in  Fig.  3.8. 
The  normalization  factor  is  WjjLchCox'  intrinsic  gate 

charge  Qqi  is  given  by  Qgj  = ~(Qgi  + Qdi  "*■  Qbi^  * Notice  that  the 
charge  quantities  are  continuous  over  all  the  device  operation 
regions.  The  differentiation  of  intrinsic  charges  with 
respect  to  the  terminal  voltages  gives  the  intrinsic  capaci- 
tance [Mas88,  Gha90,  Rho93] ; 


(3.46) 


Q / (V) 
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Fig.  3.8  Simulated  intrinsic  channel  charges  associated  with 

the  gate,  drain,  source,  and  body.  Note  that  -Qsi 
is  plotted  to  avoid  the  overlapping  of  data. 
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where  is  either  of  four  intrinsic  charges  ( Qgj,  Q^j,  Qgj, 
Qgi  ) ' is  the  voltage  of  j terminal  ( gate,  drain,  or 
source  ),  and  5^^  = -1  for  i^j  and  = 1 for  i = j . Twelve 
different  intrinsic  capacitances  can  be  defined  according  to 
(3.46) . Four  of  them  are  shown  in  Fig.  3.9.  The  capacitances 
are  continuous  at  the  boundaries  of  the  linear  and  satura- 
tion region.  In  saturation,  the  drain  voltage  is  isolated 
from  the  channel  and  cannot  affect  the  intrinsic  charges. 
Thus,  all  capacitances  are  zero.  However,  the  model  shows  an 
abrupt  change  in  the  slope  between  the  two  regions  because  of 
the  neglect  of  a slight  channel  length  modulation  effect  in 
our  current  saturation  model.  One  of  the  interesting  behav- 
iors for  the  graded -channel  charges  is  seen  in  the  Cgo  and 
Cqd  characteristics.  The  graded  channel  in  the  power  DMOST 
typically  has  a higher  doping  density  than  in  a conventional 
MOSFET.  So,  at  low  Vqs»  the  body  depletion  charge  is  more 
sensitive  to  increasing  ( i.e.,  increasing  V^h  ) than  the 
inversion  charge  ( as  shown  in  Fig.  3.8  ) . Hence  Cbd  sur- 
passes CsD  and  Cdd/  and,  accordingly,  Cqd  becomes  negative  at 
low  Vpg  in  contrast  to  conventional  MOSFETs  [Chu89,  Bud90] . 

Fig.  3.10  shows  the  characteristics  for  two  important 
device  capacitances,  the  reverse  transfer  capacitance 
( ) and  the  output  capacitance 
^^oss  “ ^ where  Vjxs  is  the  drain-to-gate  voltage,  Qqt 
^ = Qgsol~Qbi“Qsi“Qdi  ■‘■^GDOL  ^ total  gate  charge,  and 


C / (WzLchC 
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Fig.  3.9  Simulated  characteristics  of  four  intrinsic  channel 

capacitances:  Cqd,  Cqq,  Cgo/  Cg^. 


C/Wz  (fF/um) 
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5.0 


measurement 

simulation 


Vqs  “ Vss  = 0.0  V 


50.0 


Vds  (V) 


Fig.  3.10  Device  capacitance  characteristics  in  the  ND- 

LDMOST:  the  reverse  transfer  capacitance  (Cj-gg)  and 
the  output  capacitance  . 
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Qdt  ( “ Qdi  ^gdol  ^JE  ^DE  ^jc  ^Dc  ^ the  total  drain 
charge.  Agreement  is  satisfactory  over  the  entire  range  of 

Yds*  The  discrepancy  in  the  output  capacitance  at  low  V^s 
results  from  the  inaccuracy  in  the  models  for  Qjg  and  Qjc  at 
low  reverse-biased  conditions.  For  Vgg  = 0 as  in  Fig.  3.10, 
the  reverse  transfer  capacitance  comes  exclusively  from  the 
extrinsic  nonplanar  gate  overlap  capacitance  over  the  drift 
region.  Although  a constant  charge  sharing  factor  (fgp  = 0.8) 
was  used  for  calculation  of  Qqdol'  ^ good  correspondence 
between  simulation  and  measurement  is  observed. 

3 . 5 Summary 

An  analytical  and  compact  model,  accounting  parasitic 
actions  like  the  BJT  and/or  the  clamping  diode  in  the  ND- 
LDMOST  structure,  has  been  presented.  Unlike  the  conventional 
BJT  analysis  method  [Gum70] , the  parasitic  model  herein  has 
been  developed  from  the  high- injection  ambipolar  transport 
equation,  and  incorporated  into  the  whole  steady-state  device 
model.  The  model  predictions  are  well  consistent  with  the 
underlying  physics.  This  parasitic  component  model,  which  has 
the  advantage  of  directly  using  LDMOST  physical  and  process 
parameters,  enables  to  simulate  another  unique  device  charac- 
teristics for  reversely  driven  V^g  during  the  DMOST  transient 
switching. 
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A charge-based  dynamic  model  has  been  also  developed  in 
this  chapter,  which  is  represented  by  a quasi-static  model  in 
terms  of  the  integrated  charge  stored  in  the  device.  The 
expressions  are  formulated  from  the  derived  dc  models  and  used 
to  calculate  the  capacitances.  The  intrinsic  graded-channel 
charges,  developed  for  the  first  time,  siir5)ly  characterize  the 
key  features  in  the  corr^lex  graded-channel,  without  signifi- 
cant loss  of  accuracy.  is  approximated  well  by  a linear 
function  of  the  position  and  apparently  simulates  the  added 
inversion  charge  away  from  the  source  in  the  laterally  graded 
channel.  Qgi,  Qdi/  and  Qgi  clearly  show  the  dependence  of  the 
charge-voltage  characteristics  on  Vqs  for  all  possible  operat- 
ing regions  ( linear,  saturation,  cut-off,  and  accumulation 
) . In  modeling  the  extrinsic  charges  of  the  nonplanar  struc- 
ture, Qgdol  characterized  by  eirploying  a charge  sharing 
factor.  The  simulations  for  all  these  regional  charges  are 
well  consistent  with  the  theoretical  predictions,  and  the 
comparisons  between  measured  and  simulated  device  capaci- 
tances are  in  a good  agreement . 


CHAPTER  4 

MODEL  IMPLEMENTATION  IN  SPICE 
4.1  Introduction 

The  computer  package  known  as  SPICE  { Simulation  Program 
with  Integrated  Circuit  Emphasis  ) is  one  of  the  most  general 
and  powerful  circuit  analysis  programs.  Models  for  semicon- 
ductor components  such  as  diodes,  BJT's,  JFET's,  and  MOSFET's 
are  already  incorporated  into  the  existing  SPICE  package 
[Nag75,  Qua89,  HSP90] . The  lateral  DMOS  transistor,  which  is 
probably  the  most  widely  used  switching  device  in  power 
electronics,  has  not  hitherto  been  available  in  the  SPICE 
model  library.  This  chapter  presents  a practical  approach  to 
incorporating  a charge-based  semi -numerical  LDMOST  model, 
described  in  the  preceding  chapters,  into  the  commercially 
available  SPICE2G.6  program;  that  is,  generating  a modified- 
SPICE  which  allow  the  accurate  simulation  of  lateral  DMOS 
transistor  circuits.  ( From  now  on,  we  call  this  modified- 
SPICE  LADISPICE  ( LAteral  Dmos  transXstor  SPICE  ) to  distin- 
guish it  from  conventional  SPICE' s.  ) The  procedures  used  to 
incorporate  the  model  into  SPICE2G.6  and  various  techniques 
necessary  to  ensure  convergence  are  described. 
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Fig.  4.1  shows  a complete  network  representation  of  the 
physical,  semi -numerical,  and  charge-based  ND-LDMOST  model 
described  in  Chapter  2 and  3 . It  consists  of  three  linear 
resistors  ( Rbd/  ^dR'  ^sub  ) » four  static  current  sources  ( 
^DS'  ^BF'  ^BR'  ^CT  ) ' nine  dynamic  current  sources  { dQQ_ 

dQgj/dt,  dQj3j/dt,  dQgj/dt,  dQ^^QL/dt,  dQjg/dt,  dQj^/dt, 

dQ^g/dt,  dQpc/dt  ),  which  yield 

“ ^DS“^BF~^BR 

+ dQj^/dt  - dQ^^/dt  + dQjg/dt  - dQ^^/dt  (4.1) 

-dQsj/dt-dQj3i/dt  (4.2) 

“ ^BF  ^CT~  ^DS  “ 

+ dQgj/dt + dQi3E/dt-dQjg/dt  (4.3) 

^SUB  ~ ^BR  ~ ^CT  ~ (4.4) 

where  the  direction  of  terminal  currents  is  into  the  device. 
The  transient  currents  are  described  by  time  derivatives  of 
the  regional  charges.  The  charges  are  characterized  in  terms 
of  the  node  voltages  through  the  quasi-static  approximation. 
In  normal  operating  conditions,  the  body-drift  (E-B)  and 
substrate-drift  (C-B)  junctions  are  reverse-biased,  so  that 
the  voltage  drops  across  series  resistances  ( Rbd/  Rdr/ 

Rsub  ) are  negligible.  Therefore,  V^g^-Vj^g  and  V^B^Vgg-Vgs, 
accordingly  only  two  terms  dQjg/dt  and  dQj^/dt  of  the  para- 
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G 


Fig.  4.1  Complete  network  representation  of  the  ND-LDMOST 

mode 1 s . 
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sitic  BJT  model  have  an  effect  on  the  device  characteris- 
tics . 

Generally,  to  incorporate  the  device  model  into  a cir- 
cuit simulator,  the  model  must  relate  the  currents  and 
charges  as  a function  of  the  node  voltages.  Table  4.1  shows 
that  each  static  current  and  every  charge  of  LDMOST  model  is 
completely  described,  quasi-statically,  in  terms  of  the 
instantaneous  voltage  differences  V^g,  Vgg,  V^b,  and  V^b*  When 
all  of  the  elements  in  the  model  are  directly  implemented  in 
SPICE,  for  a given  Newton-Raphson  iteration  on  node  voltages 
in  SPICE  nodal  analysis,  the  model  will  provide  an  element 
current  and  the  needed  partial  derivatives  with  respect  to 
the  node  voltages . A special  note  should  be  given  on  the  dc 
current  source  I^g  and  the  intrinsic  channel  charges  ( Qgj, 
Qbj,  and  Qbi  ) , which  are  implicit  function  of  both  V^g  and 
^GS*  ios  is  based  on  integration  of  the  regional  ( channel 
and  drift  ) analysis.  As  indicated  Fig.  2.6,  for  given  Vpg 
and  Vgg,  the  model  element  numerically  solves  the  implicit 
system  of  equations  to  yield  I^s  and  V^h.  Then,  the  numeri- 
cally solved  is  used  for  calculation  of  the  intrinsic 
channel  charges  Qgj,  Q^j,  and  Qbi- 

In  order  to  implement  this  seven-node  ND-LDMOST  model 
into  SPICE2G.6  program,  new  subroutines  for  the  LDMOST  model 
have  been  added,  and  the  source  code  of  the  existing  subrou- 
tines has  been  partially  modified  to  accept  the  new  model.  In 
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TABLE  4.1 

BRANCH  VOLTAGE  DEPENDENCE  OF  MODEL  ELEMENTS 


Model  elements 

Voltage  dependence 

^BD  • ^DR ' ^SUB 

constant 

^DS 

^GS 

^BF'  ^BR'  ^CT 

^CB 

Qgsol 

VqS 

Qsi'  Qdi'  Qbi'  Qgdol 

^DS'  ^GS 

Qje 

^EB 

Qjc 

VcB 

QdE  ! Qdc 

^EB'  ^CB 
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Section  4.2,  the  overall  structure  of  software  is  described. 
The  procedures  for  subroutine  addition  and  modification  are 
described  in  Section  4.3  and  4.4,  respectively.  Section  4.5 
presents  the  input  format  of  the  new  lateral  DMOS  transistor 
model,  and  the  model  parameters  which  completely  describe  the 
static  and  dynamic  behavior  of  LDMOSTs  within  LADISPICE.  In 
Section  4.6,  the  utility  and  computing  efficiency  of  the 
LADISPICE  are  demonstrated  by  simulating  several  representa- 
tive power  semiconductor  circuits.  These  simulations  confirm 
that  the  LADISPICE  may  be  used  as  an  efficient  tool  for  the 
design  of  LDMOST  integrated  circuits. 

4.2  Overall  Structure  of  LADISPICE 

The  SPICE2G.6  program  consists  of  seven  major  indepen- 
dent modules  ( READIN,  ERRCHK,  SETUP,  DCTRAN,  DCOP,  ACAN,  and 
OVTPVT  ) together  with  a main  program  [Nag75]  . The  main 
program  SPICE. f contains  the  main  control  loop  of  the  whole 
program.  It  calls  the  various  modules  in  sequence  that  is 
required  for  the  specific  simulation.  The  READIN  module  reads 
the  input  file,  checks  the  basic  syntax  errors  in  the  input, 
and  constructs  the  data  structure  for  the  circuit . The  ERRCHK 
module  checks  common  user  errors  in  the  circuit  description. 
Then,  it  constructs  an  ordered  list  of  the  user  node  numbers, 
and  prints  the  circuit  listing,  the  device  model  parameter 
listing  and  the  node  table.  The  SETUP  module  constructs  the 
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integer  pointer  structure  that  is  used  by  both  the  DCTRAN  and 
ACAN  module.  The  largest  and  most  complicated  module  DCTRAN 
performs  the  dc  operating-point  analysis,  the  dc  transfer- 
curve  analysis,  the  transient  initial-condition  analysis,  and 
the  transient  analysis.  The  DCOP  module  prints  out  information 
about  the  device  operating  point,  computes  the  small-signal 
transfer  function,  and  computes  dc  sensitivities  of  output 
variables  with  respect  to  circuit  parameters . The  ACAN  module 
performs  the  small-signal  frequency  domain  analysis,  the 
noise  analysis,  and  the  distortion  analysis.  Finally,  the  last 
module  OVTPVT  generates  the  tabular  listings  of  the  output 
variables  and  the  line-printer  plots  of  the  simulated  results. 

Based  on  the  SPICE2G.6  structure,  three  new  subroutines 
{ LDMOST.f,  LDTMOS.f,  and  LDTBJT.f  ) have  been  added  for  dc 
and  transient  analysis  of  the  lateral  DMOS  transistor,  and  the 
source  codes  of  18  original  subroutines  in  SPICE2G.6  have  been 
modified  appropriately  to  accept  the  new  model,  as  listed  in 
Table  4.2.  { Hence,  LADISPICE  is  coir^osed  of  total  127  subpro- 
grams: one  main  program  SPICE. f and  126  subroutines.  ) The 
model  for  lateral  DMOS  transistor  is  accessed  similarly  to 
other  semiconductor  device  models  written  into  SPICE2G.6.  The 
LDMOST  is  listed  with  an  identifying  number,  say  ID  = 15,  of 
the  SPICE  list  contents  with  the  letter  Y as  its  prefix,  and 
the  LDMOST  device  model  is  listed  with  ID  = 25.  ( The  detailed 


87 


TABLE  4.2 

NEW/MODIFIED  SUBROUTINES  IN  SPICE2G.6  SOURCE  CODE 

TO  IMPLEMENT  LDMOST  MODEL 


Subroutine  name 

New  subroutines 

LDMOST . f , 

LDTMOS.f,  LDTBJT.f 

Modified  subroutines 

ACLOAD . f , 

ADDELT . f , 

ALTER . f , 

DCOP.f,  DCTRAN.f, 

ELPRNT . f , 

ERRCHK.f,  FIND.f, 

LNKREF . f , 

LOAD.f,  MATLOC.f, 

MATPTR . f , 

MODCHK . f , 

OUTDEF . f , 

READIN . f , 

SPICE. f. 

TOPCHK.f,  TRUNC.f 
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linked-list  specifications  for  the  lateral  DMOS  transistor 
are  given  in  Appendix  A.  ) 

Fig.  4.2  shows  the  calling  diagram  of  the  LDMOST  model 
routine  for  dc  and  transient  analysis,  with  other  semiconduc- 
tor device  routines  ( DIODE. f,  BJT.f,  JFET.f,  andMOSFET.f  ). 
For  the  specific  analysis  requested  by  the  user,  the  subrou- 
tines DCTRAN.f  and  LOAD.f  call  the  subroutine  LDMOST. f,  and 
subsequently  LDMOST. f calls  the  separate  routines  LDTMOS.f 
and  LDTBJT.f,  within  the  module  DCTRAN.  The  subroutine  ITERS. f 
determines  the  solution  of  the  circuit  equations.  The  subrou- 
tine LOAD.f  constructs  the  linearized  system  of  nodal 
equations  for  the  specific  iteration.  The  contributions  of 
admittance  matrix  from  the  LDMOSTs  are  computed  separately  in 
the  subroutine  LDMOST. f.  ( The  detailed  nodal  admittance 
matrix  for  the  LDMOST  is  shown  in  Appendix  B.  ) 

Basically,  the  LDMOST. f routine  links  the  LDMOST  model 
to  the  SPICE  nodal  analysis.  The  structure  of  LDMOST. f is 
similar  to  that  of  the  original  semiconductor  device  routines 
in  SPICE2G.6  except  that  most  of  the  model  calculations  are 
done  in  the  separate  routines  LDTMOS.f  and  LDTBJT.f.  The 
LDTMOS.f  routine  calculates  the  model  elements  for  the  normal 
operating  conditions:  I^g,  Qsi/  Qdi»  Qbi»  Qgdol'  returns 
their  calculated  values  to  the  LDMOST. f routine.  Similarly, 
the  LDTBJT.f  routine  evaluates  the  currents  ( Ibf»  ^br'  ^ct  ) 
and  charges  ( Opg,  Qj^;  ) related  to  the  parasitic  BJT  action. 
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Fig.  4.2 


Calling  diagram  of  lateral  DMOS  transistor  model 
routine  with  the  other  semiconductor  device  rou- 
tines . 
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and  returns  the  calculated  values  into  the  subroutine 
LDMOST.f.  The  other  terms  ( Qqsol'  Qje'  Qjc  ) model  are 
calculated  in  the  inside  of  the  subroutine  LDMOST.f  for  the 
nvimerical  efficiency.  Accordingly,  the  LDMOST  model  equations 
are  completely  calculated  and  linked  to  the  SPICE  nodal  analy- 
sis by  the  new  added  three  subroutines. 

4.3  Addition  of  New  Subroutines 

This  section  presents  the  detailed  structure  of  new  added 
model  routines  ( LDMOST.f,  LDTMOS.f,  LDTBJT.f  ) using  the 
flowchart.  { The  complete  source  codes  for  these  3 subroutines 
are  given  in  Appendix  C.  ) The  numerical  techniques  used  to 
solve  the  model  equations  are  discussed,  and  the  redef inements 
of  model  equations  to  resolve  the  convergence  problem  in  SPICE 
are  described. 

4.3.1  Subroutine  LDMOST . f 

As  the  first-level  device  routine  for  dc  and  transient 
simulation  of  lateral  DMOS  transistors,  the  LDMOST.f  routine, 
shown  in  Fig.  4.3,  connects  the  model  to  the  SPICE  nodal 
analysis.  LDMOST.f,  first,  initializes  local  variables  and 
model  parameters  used  in  the  model  routine.  The  extrinsic 
resistances  Rbd»  ^dr»  ^sub»  which  are  preprocessed  in 
MODCHK.f,  are  treated  as  conductance  elements  in  LDMOST.f; 
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Fig.  4.3  Flowchart  of  the  subroutine  LDMOST.f. 
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that  is,  their  reciprocal  values  are  directly  added  into  the 
admittance  matrix.  The  node  voltage  differences  V^g,  VQg,  Vgg, 
and  VcB  are  initialized  and  computed  based  on  the  user's 
initial  conditions  read  through  input  file  or  their  previous 
values.  Then,  those  branch  voltages  are  limited  by  calling 
subroutines  FETLIM.f,  LIMVDS.f,  and  PNJLIM.f  to  ease  the 
convergence . 

Next,  LDMOST.f  calls  successively  LDTMOS.f  and  LDTBJT.f 
which  calculate  the  operating-point  currents  and  charges,  for 
the  normal  mode  and  the  parasitic  BJT  mode,  respectively.  For 
each  call,  the  two  internal  node  voltage  differences  ( Vgg  and 
Vq3  for  normal  mode,  Vge  and  Vcb  for  parasitic  BJT  mode  ) are 
passed  into  LDTMOS.f  and  LDTBJT.f,  and  the  needed  current 
components  and  charge  coirponents  are  calculated  and  returned 
to  LDMOST.f.  Because  of  implicit  characteristics  of  model 
equations,  the  partial  derivatives  of  currents  and  charges 
needed  in  the  SPICE  nodal  analysis  are  evaluated  by  these 
multiple  calls  of  the  LDTMOS.f  and  LDTBJT.f  routines  with 
perturbed  values  of  branch  voltages,  V^g,  Vgs»  Veb»  and  Vqb* 
These  derivatives,  that  is,  conductances  and  capacitances  are 
calculated  by  the  second-order  (forward-difference)  numerical 
differentiation. 

Then,  LDMOST.f  performs  the  jobs  for  evaluating  the 
equivalent  currents  needed  in  the  admittance  matrix.  For  the 
dc  analysis,  the  equivalent  charge  currents  are  simply  set  to 
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zero.  During  the  transient  analysis,  LDMOST.f  calculates  the 
explicit  charges  ( Qqsol'  Qje'  Qjc  ) their  capacitances, 
respectively.  Then,  LDMOST.f  sequencially  proceeds  evaluating 
the  equivalent  conductances  for  all  of  the  capacitances,  and 
calculating  the  equivalent  charge  currents. 

Lastly,  LDMOST.f  calculates  the  equivalent  dc  currents 
based  on  conductances,  loads  the  elements  of  circuit  admit- 
tance matrix,  and  then  returns  the  control  to  the  SPICE  nodal 
analysis . 

4.3.2  Subroutine  LDTMOS . f 

The  task  of  the  LDTMOS. f routine  is  to  return  the  drain- 
source  current  I^s,  the  intrinsic  graded-channel  charges  ( 
Qsi'  Qdi»  Qbi  ) ' the  gate-drain  overlap  charge  Qqdol  the 
two  input  branch  voltages  Vog  and  Vgg.  The  resulting  structure 
of  LDTMOS. f is  shown  in  Fig.  4.4.  First,  the  local  variables 
and  model  parameters  are  initialized  based  on  the  pointer  LOG 
passed  from  the  LDMOST.f  routine.  All  the  physical  quantities 
that  can  be  calculated  outside  the  iterative  loop  are  done  for 
the  numerical  efficiency.  Then,  the  I^g  model  solution  loop 
begins,  which  is  same  as  the  semi -numerical  routine  in  Fig. 
2.6.  For  the  normal  operating  bias  conditions  ( Vpg>0  and 
Vq2>Vt  ),  an  initial  guess  is  made  for  the  drift  region 
voltage  drop  and  the  channel  current  is  calculated  and 
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Flowchart  of  the  subroutine  LDTMOS.f.  and  V*Qf^ 
are  the  guessed  value  and  computed  value,  respec- 
tively, for  drift  region  voltage  drop. 


Fig.  4.4 
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set  to  be  equal  to  the  drain-source  current  either  the 

linear  or  saturation  region. 

The  next  step  is  the  drift  region  analysis.  A system  of 
differential  equations  in  the  equations  (2.28)  and  (2.29) 


are  solved  by  the  Runge-Kutta  method  for  one  of  11  possible 
2-D  current  paths,  which  is  determined  by  the  two  terminal 
biases  V^s  and  Vqs.  The  R-K  method  is  numerically  explicit, 
and  evaluates  the  functions  at  selected  points  on  each 
subinterval.  In  this  worlc,  the  effective  1-D  length  of  current 
path  WgcR  is  divided  by  a constant  n4,  so  that  the  increment 
h^  = Wg(,j^/n4 . For  given  initial  conditions  for  potential, 
electric  field,  and  coordinate;  ~ ^ch ' ~ 

u(®)  =0,  respectively,  the  solutions  at  the  next  point  are 
obtained  by  the  following  algorithm  [Con80] : 


g4i(u,E^)  = 


ho  + utany  ®J  Eg 


1 ^0  1 . 
Dn  H 


(4.5) 


942 


(4.6) 


u (n)  = n • h^ 


(4.7) 


)Cj  +2)^2  + 2)C3  + 1^4 


(4.8) 


6 


y(n  + 1)  = y(n) 


IJ  + 2I2  + 2I3  + I4 


(4.9) 


6 


= ^4-  g^Ju<"),E^'^)] 


(4.10) 
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ll  = h,-g,2[E^(n)] 

^2  = N • 941  + T’  y] 

^2  = W'942[E^"^+y] 

^3  =~^4‘  941  y ’ + y] 

13  =N-942[E,(")+y] 

^4  = h4  • [u<")  + h4,  E^n)  + kj] 

14  = h4'942[E^(^>+k3] 

where  n = 0,  1,  2,  ...,  n4-l.  Therefore,  is  obtained  by 

the  potential  difference  between  u = Wg^,ji  and  u = 0:  that  is, 
V*pfi  = . This  Runge-Kutta  method  tends  to  increase 

the  simulation  time,  but  gives  great  accuracy  and  avoids 
solving  the  complex  system  of  differential  equations 
analytically . 

If  the  calculated  V*qr  is  not  close  to  the  guessed  V^p. 
within  a specified  tolerance,  the  new  guess  for  Vqr  is 
estimated  using  the  Secant  method  [McC88] : 


(4.11) 

(4.12) 

(4.13) 

(4.14) 

(4.15) 

(4.16) 

(4.17) 


(j+i) 

DR 


(4.18) 


(j) 

DR 


(4.19) 
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where  j is  an  integer.  Otherwise,  the  Ips  model  solution  loop 
is  exited.  Then,  the  regional  charges  Qgj,  Qqj,  Qgj,  and  Qqdol 
are  calculated  in  sequence. 

Generally,  in  order  to  resolve  the  convergence  problems 
in  SPICE,  the  model  equations  should  be  continuous  for  all  of 
the  terminal  bias  variations.  Several  attempts  were  already 
made  to  eliminate  any  discontinuity  by  redefining  pertinent 
model  equations.  In  Eq.  (2.9),  the  VQg-dependent  kf3  was  intro- 
duced to  prevent  a discontinuity  of  Iqh  at  V^g  = V,p . In  Eq. 
(2.20),  yfox  ^as  redefined  to  ensure  continuous  vertical 
movement  of  line-4.  Another  excimple  is  setting  I^g  = 0 for 
Vog<0,  described  in  Fig.  2.6. 

In  addition  to  these  redefinements,  a few  more  attempts 
are  made  in  LDTMOS . f for  the  robust  model  for  all  the  bias 
regions.  During  the  iteration  of  the  I^g  solution  loop  in  Fig. 
4.4,  Vqr  can  be  guessed  to  a value  greater  than  V^g,  i.e.,  Vch 
can  be  negative.  In  this  case,  a conputational  logic  error  may 
appear  in  some  model  equations.  To  prevent  this  error,  artifi- 
cial model  equations  for  are  added  in  the  mobility  term 
at  Eq.  (2.3),  and  in  the  boundary  condition  for  the  electric 
field  at  Eq.  (2.31): 


= >‘„ch(>-Vch) 


(for  V .SO)  (4.20) 


(for  V.^SO)  (4.21) 
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Additionally,  for  the  bias  conditions  of  VQg>V,p  and  V^g^O, 
the  intrinsic  channel  charges  ( Qsi»  Qdi'  Qbi  ) set  to  the 

values  of  V = 0 in  their  linear  region  model  equations 
(3.24),  (3.25),  and  (3.29).  This  removal  of  model  discontinu- 
ities and  the  complete  model  equations  covering  all  of  the 
bias  regions  are  quite  important  for  convergence  in  SPICE. 

4.3.3  Subrout ine  LDTB JT . f 

The  purpose  of  the  LDTBJT.f  routine  is  to  evaluate  the 
currents  ( Igp,  Igp,  ) and  charges  ( Q^g,  Qqc  ) related  to 
the  parasitic  BJT  action  for  the  two  input  branch  voltages  Vge 
and  VcB/  and  to  return  the  calculated  values  into  the  subrou- 
tine LDMOST.f.  Because  the  parasitic  BJT  model  equations  are 
fully  analytical  with  respect  to  the  internal  node  voltage 
differences  Vge  and  Vcb»  the  structure  of  LDTBJT.f,  shown  in 
Fig.  4.5,  is  similar  to  that  of  LDTMOS.f  except  that  LDTBJT.f 
doesn't  need  the  internal  iteration  loop. 

A feature  in  the  subroutine  LDTBJT.f  is  to  avoid  the 
numerical  underflow  in  LADISPICE  simulation.  The  underflow 
problem  can  occur  in  all  exponential  terms  involving  negative 
voltages  in  their  exponents.  For  exaitple,  in  the  normal 

operating  conditions  ( '^eb-~^ds  ^cb-^ss~^ds 

exponential  terms  involving  V^b  and  Vcb  should  be  zero  for 
large  Vqs*  However,  they  are  numerically  nonzero,  being  a very 
small  mamber  which  may  not  be  defined  in  a given  precision  of 
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Fig.  4.5 


Flowchart  of  the  subroutine  LDTBJT.f. 
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a specific  computer  system.  This  underflow  problem  is  resolved 
by  limiting  the  negative  exponents  to  some  value  which  does 
not  affect  the  model  characteristics.  In  LDTBJT.f,  this  limit- 
ing is  done  by  forcing  the  exponential  terms  to  zero  if  Vgg 
and  Vqb  in  their  exponents  are  less  than  -200(kT/q)  . 

As  in  LDTMOS.f,  some  redef inements  of  model  equations  are 
also  done  in  LDTBJT.f  to  ensure  the  model  continuity  for  all 
the  bias  regions.  The  parasitic  BJT  model  described  in  Chapter 
3 is  valid  only  for  Wgjg>0;  that  is,  it  does  not  characterize 
the  punch-through  effect  in  the  base.  Although,  we  don't  need 
to  consider  it  for  practical  bias  conditions  in  the  long  base- 
width  device  like  parasitic  BJT  of  the  LDMOST,  artifitial 
model  equations  are  added  in  LDTBJT.f  for  the  robustness  of 
the  model.  First,  the  term  related  with  recombination  in  the 
quasi-neutral  base  region 

cosh  (W,^)  - 1 

sinhCWj^)  ' ' • ' 

appearing  in  Igp,  Igp,  Qqe,  and  Qix:  ( Eqs.  (3.8),  (3.9),  (3.42), 
and  (3.43)  ),  is  forced  to  zero  for  Wgjj^O,  which  does  not 

hurt  the  BJT  physics.  The  other  redef inement  is  done  for  the 
term  expressing  the  Early  effect  in 

1 + cosh  (Wg^)  /bg 

sinh(Wgjj) 


(4.23) 
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which  increases  to  infinite  as  Wq^  approaches  to  zero,  imply- 
ing the  punch- through.  It  is  not  defined  for  Wgjj^O.  This 
problem  is  resolved  by  limiting  the  term  at  a constant  value 
with  Wgjj  = 10~^  if  Wbn  is  less  than  10“^;  that  is,  the  quasi- 
neutral base  width  is  less  than  0.1  % of  the  base  ambipolar 
diffusion  length. 

4.4  Modification  of  Original  Subroutines 

With  addition  of  the  new  model  subroutines,  the  source 
code  in  the  original  SPICE2G.6  program  has  been  partially 
modified.  Following  the  detailed  procedures  in  [Fit 91] , 18 
existing  subroutines  ( see  Table  4.2  ) have  been  changed  in 
order  to  accept  the  developed  LDMOST  model . 

Subroutine  SPICE. f was  modified  to  change  the  print 
header  information  reflecting  the  new  version  of  SPICE,  and 
to  include  LDMOSTs  in  the  job  statistics  summary.  Subroutine 
READIN.f,  which  deals  with  input  information,  was  modified  in 
order  that  the  proposed  new  format  of  the  LDMOST  can  be  read 
and  any  missing  parameter  can  be  checked.  In  subroutine 
MODCHK.f,  the  source  code  was  changed  to  perform  a preprocess- 
ing of  LDMOST  model  parameters,  to  print  out  the  device  model 
summary,  and  to  reserve  additional  internal  nodes  required  by 
extrinsic  resistances  Rbd»  ^dr'  ^sub*  *^he  source  code  of 
subroutine  ERRCHK.f  was  also  changed  to  allow  that  the  node 


initial  conditions  of  LDMOSTs  can  be  translated  into  device 
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initial  conditions,  and  that  the  default  values  for  LDMOST 
geometries  can  be  used  if  the  user  does  not  specify. 

Subroutines  ADDELT.f,  MATLOC.f,  and  MATPTR.f  were 
modified  to  add  LDMOSTs  in  the  circuit  definition  lists,  and 
to  set  up  and  store  the  nonzero  coefficient  matrix  for  the 
LDMOSTs.  Subroutine  LOAD.f  was  also  modified  to  process 
LDMOSTs  in  the  loading  operation  of  the  admittance  matrix. 
These  modified  subroutines  then  build  up  the  system  matrix 
representing  the  simulated  circuit. 

In  subroutine  DCTRAN.f,  a calling  statement  for  LDMOSTs 
was  added  and  a local  variable  (NUMESE)  was  expanded,  to 
ensure  that  the  dc  and  transient  analysis  are  correctly 
performed  when  LDMOSTs  are  in  the  circuit.  In  subroutines 
ELPRNT.f  and  DCOP.f,  sections  of  code  were  added  to  print  out 
the  user  node  number,  the  device  width,  and  the  geometric 
areas  ( Ag,  A^,  Ag^  ) for  LDMOST  model,  as  well  as  to  allow 
printing  out  the  operating  point  variables  for  LDMOSTs. 

Other  existing  subroutines  such  as  ACLOAD.f,  ALTER. f, 
FIND.f,  LNKREF.f,  TOPCHK.f,  TRUNC.f,  and  OUTDEF.f  were  appro- 
priately modified  in  order  that  the  LADISPICE  might  become  a 
general-purpose  simulator  for  electronic  circuits. 

4.5  Input  Format  and  Model  Parameters 

In  LADISPICE,  the  input  format  of  the  LDMOST  model  is  the 
free-format  type  and  essentially  similar  to  that  of  the  other 
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semiconductor  devices  in  SPICE2G.6,  shovm  in  Fig.  4.6.  The 
device  card  starts  with  a Y-prefix  element  name.  ND,  NG,  NS, 
and  NSUB  are  the  drain,  gate,  source,  and  substrate  node, 
respectively,  and  MNANE  is  the  model  name.  The  layout-depen- 
dent  parameters  WZ,  AE,  AC,  and  AEC,  which  can  differ  from 
device  to  device  on  the  same  power  IC  chip,  are  set  by  the 
user  { the  definitions  and  default  values  are  tabulated  in 
Table  4.3  ) . OFF  indicates  an  optional  initial  condition  on 
the  device  for  dc  analysis.  At  the  end  of  the  device  card, 
initial  conditions  for  device  terminal  voltages  can  be 
supplied  as  another  option.  These  initial  conditions  can  be 
useful  in  easing  convergence  difficulties.  The  .MODEL  card 
contains  the  model  name  and  the  device  type  NLDMT,  indicating 
the  n-channel  lateral  DMOS  transistor.  Subsequently,  it 
contains  a set  of  36  model  line  parameters,  of  which  defini- 
tions and  default  values  are  listed  in  Table  4.4  { The 
definition  of  device  geometric  parameters  is  indicated  in  Fig. 
1.1.  ).  These  parameters  are  in  general  common  to  all  the 
LDMOST  devices  in  the  same  power  IC  chip,  and  are  closely 
related  to  process  information  and  the  device  physics. 

Basically,  simulation  accuracy  is  critically  dependent 
on  the  proper  choice  of  model  parameters . These  parameters  can 
be  fitted  to  a particular  fabrication  process  through  exper- 
imental characterization  of  the  test  element  group  (TEG) . So, 
effective  use  of  the  circuit  simulator  demands  a well-devel- 
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YxxxxxUD  NG  NS  NSUB  MNAME  <}NZ=vat>  <AE=val> 
<AC^val>  <AEC—val>  <OFF>  <IC=  Vgg> 


.MODEL  MNAME  NLDMT  ( para1=va/,  para2=va/, ... ) 


Fig.  4.6  Input  format  of  the  lateral  DMOS  transistor  model 

in  LADISPICE. 


TABLE  4.3 

LATERAL  DMOS  TRANSISTOR  DEVICE  LINE  PARAMETERS 


Name 

Description 

Units 

Default 

WZ 

Device  width 

m 

l.Oe-6 

AE 

Effective  body-drift 
junction  area 

rr.2 

m 

5.0e-12 

AC 

Effective  substrate-drift 
junction  area 

m2 

m 

l.Oe-11 

AEC 

Effective  parasitic  BJT 
cross-sectional  area 

m2 

m 

• 

O 

(D 

1 

to 
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TABLE  4.4 

LATERAL  DMOS  TRANSISTOR  MODEL  LINE  PARAMETERS 


Name 

Description 

Units 

Default 

LCH 

Channel  length 

cm 

l.Oe-4 

LDl 

Overlapped  length  of  thin  gate 
oxide  over  drift  region 

cm 

1 . Oe-4 

LD2 

Tapered  oxide  length 

cm 

0.6e-4 

LD3 

Overlapped  poly-gate  length 
over  thick  field  oxide 

cm 

2.5e-4 

DELTAL 

Lateral  source  diffusion 
length 

cm 

0.35e-4 

TOX 

Thin  gate  oxide  thickness 

cm 

5 . Oe-6 

TFOX 

Thick  field  oxide  thickness 

cm 

1 

0) 

CO 

• 

o 

YJB 

Body  junction  depth 

cm 

2 . Oe-4 

YJDP 

Deep-p'*'  body  junction  depth 

cm 

3 .5e-4 

YEPI 

Drift-region  epitaxial 
thickness 

cm 

• 

O 

CD 

1 

NAO 

Peak  body-doping  density  at 
source  end 

cm"^ 

5.0el7 

NDD 

Drift-region  doping  density 

cm“^ 

l.Oeie 

NAS 

Substrate  doping  density 

-3 

cm 

5.0el5 

NABEFF 

Effective  deep-p"^  body-doping 
density 

-3 

cm 

7.0el7 

VSAT 

Electron  saturation  velocity 

cm/ sec 

1.0e7 

THETA 

Channel  mobility  degradation 
factor  due  to  vertical 
electric  field 

v-i 

1.5e-l 

UNCHO 

Low-field  electron  mobility 
in  channel 

cxvr  /V- sec 

6 . 0e2 

UND 

Low-field  electron  mobility 
in  drift  region 

ever /V- sec 

1.3e3 

BETAN 

Channel  mobility  degradation 
factor  due  to  lateral  electric 
field 

to 

• 

o 
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TABLE  4.4  --  continued 


Name 

Description 

Units 

Default 

BB 

Electron-to-hole  mobility 
ratio  in  drift  region 

— 

2.0 

TAUNE 

Electron  lifetime  in  deep-p"^ 
body 

sec 

O.le-6 

TAUHB 

High  injection  carrier 
lifetime  in  drift  region 

sec 

l.Oe-6 

TAUHC 

High  injection  carrier 
lifetime  in  substrate 

sec 

l.Oe-6 

LAB 

Ambipolar  diffusion  length  in 
drift  region 

cm 

3 . Oe-3 

LAC 

Ambipolar  diffusion  length  in 
substrate 

cm 

3 .Oe-3 

VFBCH 

Flatband  voltage  in  channel 

V 

-1.0 

VFBDR 

Flatband  voltage  in  drift 
region 

V 

-0.5 

PSIBD 

Built-in  potential  of 
body-drift  junction 

V 

0.7 

FQD 

Fractional  parameter  of 
gate-drift  overlapped  charge 

— 

0.8 

KFO 

Coefficient  of  average  body 
depletion  capacitance 

— 

0.1 

KFl 

Coefficient  of  initial 
channel  thickness  in  drift 
region 

3.0 

KF2 

Coefficient  of  lateral 
electric  field  at  channel  end 

— 

1.0 

KF4 

Coefficient  in  onset  voltage 
for  channel  current  saturation 

— 

0.9 

RBD 

Deep-p'^  ohmic  resistance 

ohm 

10.0 

RDR 

Drift-region  ohmic  resistance 

ohm 

50.0 

RSUB 

Substrate  ohmic  resistance 

ohm 

100.0 
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oped  parameter  extraction  scheme.  There  are  several 
approaches  to  parameter  extraction  for  a SPICE-based  device 
model.  Generally,  model  parameters  are  categorized  into  small 
groups  based  on  the  nature  of  parameters,  and  each  group  of 
parameters  is  extracted  sequentially  assuming  the  rest  of  the 
parameters  are  accurate.  Then,  in  order  to  produce  the  best 
fit  of  the  model  to  the  measured  data  which  are  in  the  operat- 
ing regions  of  interest,  all  of  parameters  are  optimized 
simultaneously  by  numerical  method.  Here,  a systematic  way  of 
determining  the  dc  and  transient  model  parameters  from  device 
information,  I-V  and  C-V  measurements  is  outlined  for  the 
lateral  DMOS  transistor. 

In  Table  4.5,  a corrqplete  set  of  40  model  parameters 
including  4 device  line  parameters  is  categorized  into  7 
groups  based  on  the  nature  of  parameter  extraction.  First,  the 
geometry  and  dopant  dependent  parameters  in  Group-1  are 
obtained  from  the  layout  and  process  information.  The  non- 
actual doping  density  NABEFF  ( » 8.0x10^”  cm"^  [Fos81]  ) and  the 
noncritical  electrical  parameter  PSIBD  ( »0.7  V ) can  be 
assumed  respectively.  Then,  AEG  and  KFl  are  estimated  from 
PISCES  simulations.  Next,  the  remaining  parameters  are 
extracted  to  fit  the  measurements.  The  Vrp-related  parameter 
VFBCH  is  extracted  to  fit  Id~Vqs  data  at  a low  V^g.  FQD  and 
VFBDR  are  extracted  to  fit  the  measurement  of  the  reverse- 
transfer  capacitance  characteristics.  Then,  in  Group-6,  five 
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TABLE  4.5 

CATEGORIZATION  OF  LDMOST  MODEL  PARAMETERS 

BASED  ON  EXTRACTION  NATURE 


Group 

Model  parameters 

Needed  information 
for  extraction 

1 

WZ,  LCH,  LDl,  LD2,  LD3 , 
DELTAL,  TOX,  TFOX, 

YJB,  YJDP,  YEPI,  AE, 
AC,  NAO,  NDD,  NAS 

layout  and  process 

2 

NABEFF,  PSIBD 

— 

3 

AEC,  KFl 

PISCES  simulation 

4 

VFBCH 

Id-Vgs  curve  at  low  Vqs 

5 

FQD,  VFBDR 

Crss'^DS  curve  with 
Vgs=Vss=0 

6 

VSAT,  THETA,  UNCHO , 
UND,  BETAN,  KFO,  KF2, 
KF4 

Id~Vds  curve  for 
different  Vgs 

7 

BB,  TAUNE,  TAUHB, 
TAUHC,  LAB,  LAC,  RBD, 
RDR,  RSUB 

IgUB  curves 

for  negative  with 

different  Vgs 
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critical  mobility  parameters  and  three  electrical  parameters 
( KFO,  KF2,  KF4  ) are  extracted  to  fit  the  measurement  of  Iq- 
Vds  characteristics  with  varying  Vqs*  Finally,  in  Group-7, 
three  extrinsic  resistances  RBD,  RDR,  RSUB  ( otherwise,  can 
be  estimated  from  process /layout  information  ) and  six  diffu- 
sion parameters  are  extracted  to  fit  the  measurement  of 
parasitic  BJT  characteristics;  that  is,  IguB  and  Ip  curves  for 
negative  V^g  with  varying  Vgg.  This  extraction  scheme  might  be 
iir^lemented  in  a general-purpose  parameter  extraction/optimi- 
zation program  ( UTMOST,  TECAP,  or  IC-CAP  ) , and  could  be 
exclusively  used  for  extraction  and  optimization  of  the  new 
LDMOST  model  parameters . 

4.6  Simulation  Capability  of  LADISPICE 

The  utility  of  a circuit  simulation  program  is  ultimately 
determined  by  applying  the  program  to  practical  circuits. 
Here,  the  final  version  of  the  simulator,  LADISPICE-1 . 1,  is 
demonstrated  for  transient  LDMOST  circuit  simulations.  These 
simulations  will  exemplify  well  how  the  LADISPICE  can  be  used 
effectively  as  a reliable  tool  for  advanced  power  IC  TCAD.  The 
circuit  chosen  for  simulations  is  a resistor-inductor  load 
switching  circuit  with  resistive  gate  drive  [Hef91] . In 
addition,  more  complex  circuit,  a single-phase  full-bridge 
DC-to-AC  converter  [Ras88,  Moh89] , is  also  chosen  for 
LADISPICE  simulation.  In  all  the  following  simulation 
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examples  of  this  section,  the  lateral  DMOS  transistors  have 
been  chosen  to  have  the  same  layout  dimension:  that  is, 

= 5x10"*  p.m,  Ag  = 0.25x10"^  m^,  A^,  = 0.55x10"^  m^,  Ag^,  = 0.2x10  ^ 
m^,  and  the  model  line  parameters  are  same  as  those  discussed 
in  Chapters  2 and  3 . 

Fig.  4.7  shows  the  circuit  diagram  of  the  inductive- 
resistive  switching  circuit  and  the  LADISPICE  input  file  for 
transient  simulation.  Note  that  in  the  .OPTIONS  card,  we  do 
not  use  any  user  option  easing  convergence  difficulties  ( 
RELTOL,  ABSTOL,  VNTOL,  TRTOL,  ITLl,  ITL4,  ITL5,  etc.  ).  They 
are  all  set  to  their  default  values.  ( This  situation  is  same 
in  the  other  examples  of  this  section.  ) Fig.  4.8  shows  the 
simulated  gate  voltage,  drain  voltage,  gate  current,  and  drain 
current  for  the  switching  circuit  in  Fig.  4.7.  The  load  resis- 
tance of  20  Q results  in  a steady-state  current  of  2.5  A for 
the  supply  voltage  of  50  V.  All  the  switching  waveforms  are 
in  accord  with  the  theoretical  prediction. 

For  the  same  circuit  of  Fig.  4.7,  Figs.  4.9  and  4.10  show 
the  simulated  gate-drive  characteristics  and  turn-off  charac- 
teristics, with  different  values  of  gate  resistance  and  load 
inductance,  respectively.  In  Fig.  4.9,  the  simulations  are 
made  for  gate  drive  resistances  that  range  from  0.2  kfi 
through  1.0  kfi  . The  gate  pulse  cimplitude  is  14  V with  50  nsec 
rise/fall  time.  The  plateaus  on  the  rising/falling  edges  in 
Vqs  waveform  is  caused  by  the  feedback  capacitance  in  the  DMOS 
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LDMOST  SWITCHING  CIRCUIT  WITH  RESISTOR-INDUCTOR  LOAD  AND  RESISTIVE  GATE  DRIVE 

Y1  2 10  0 NMODE  WZ=0.05  AE^.25U  ACsO.SSU  AEC=0.2U 

RG  1 3 0.5K 
RL  2 4 20 
LL  4 5 0.5U 

VDD  5 0 DC  50 

VIN  3 0 PULSE  (0  14  0 SON  SON  400N  900N ) 

.MODEL  NMODE  NLDMT 

-I-  (LCH=0.84E-4  LD1-1.0E-4  LD2=0.4E-4  LD3=2.6E-4  DELTALs0.35E-4 
+ TOX^7.0E-6  TFOXs0.5E-4  YJB=1.7E-4  YJDP=3.2E-4  YEPI=7.2E-4 
+ NA0=3.0E17  NDD=1.5E16  NAS=5.0E15  NABEFF=8.0E17  VSAT=1.0E7 
+ THETA=1.97E-1  UNCH0=4.0E2  UND=1.2E3  BETAN=1.S  BB=3.5 
+ TAUNE=0.15E-6  TAUHB=0.3E-6  TAUHC=0.3E-6  LAB=8.6E-4  LAC=23.6E-4 
+ VFBCH=-2.57  VFBDR=-0.2  PSIBD=0.7  FQD=0.8  KF0=0.1  KF1=3.0 
+ KF2=1.0  KF4=0.94  RBD=3.5  RDR=56.7  RSUB=62.1  ) 

.TRAN  ION  2100N  300N 

.PRINT  TRAN  V(3)  V(1)  V(2)  l(VDD)  l(VIN) 

•PLOT  TRAN  V(3)  V(1)  V(2)  l(VDD)  l(VIN) 

.OPTIONS  LIST  NODE  ACCT  LIMPTS=1001 

.WIDTH  IN=80  OUT=80 

.END 


Fig.  4.7  Circuit  diagram  of  resistor-inductor  load  switching 

circuit  with  resistive  gate  drive  and  input  file  for 
LADISPICE  simulation. 
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Fig.  4.8  LADISPICE-simulated  switching  waveforms  of  resis- 
tor-inductor load  switching  circuit  with  resistive 
gate  drive:  (a)  voltages  and  (b)  currents. 
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Fig.  4.9  LADISPICE-simulated  gate  characteristics  of  resis- 
tor-inductor load  switching  circuit  with  resistive 
gate  drive:  (a)  voltages  and  (b)  currents. 
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Fig.  4.10  LADISPICE-simulated  turn-off  characteristics  of 

resistor-inductor  load  switching  circuit  with 
resistive  gate  drive. 
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transistor  [She93] . The  predicted  trends  for  different  gate 
resistance  values  are  in  qualitatively  accord  with  the  results 
in  references  [Hef91,  She93,  Met94]  . This  agreement  validates 
our  nonplanar  gate  capacitance  model  in  the  investigated  DMOS 
structure.  In  Fig.  4.10,  the  turn-off  characteristics  of  drain 
voltage  and  drain  current  are  simulated  with  the  small  gate 
resistance  { 50  ) . This  is  in  order  to  exclude  the  plateaus 
effect  in  turn-off  characteristics.  The  load  inductances  are 
chosen  so  that  the  drain  voltage  overshoot  will  not  exceed 
twice  of  the  supply  voltage  ( Vdd  = 50  V ) . The  gate  voltage 
pulse  amplitude  of  14  V is  same  as  that  in  the  simulation  of 
the  gate-drive  characteristics  and  the  rise/fall  time  of  the 
gate  pulse  generator  is  also  50  nsec.  The  predictions  of  the 
current/voltage  turn-off  characteristics  with  different 
inductor  values  adequately  describe  the  overshoot  for  a given 
inductance  [Hef90,  She93] . 

Now  the  use  of  LADISPICE  is  exemplified  in  the  device 
design  for  the  switching  circuit,  based  on  the  typical  lateral 
DMOS  technology.  The  goal  of  this  example  is  to  show  the  trend 
in  circuit  performance  depending  on  process  technology.  The 
LADISPICE-predicted  turn-on  characteristics  ( the  definitions 
of  switching  times  are  inside  of  the  figure  ) versus  body- 
diffusion  junction  depth  is  plotted  in  Fig.  4.11,  with  assump- 
tion that  the  lateral  diffusion  length  of  source/body  is 
approximately  70  % of  their  junction  depths.  The  simulations 
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Fig.  4.11  LADISPICE-simulated  turn-on  switching  time  of 

resistive  load  switching  circuit.  The  input  gate 
pulse  amplitude  is  10  V with  20  nsec  rise/fall 
time . 
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are  made  for  circuit  parameters  similar  to  those  used  in  the 
reference  [Har92]  ( Rj^  = 30  Q,  = 10  Q and  the  load  induc- 
tance is  set  to  zero  ) . The  supply  voltage  is  50  V,  and  the 
gate  pulse  amplitude  is  10  V with  20  nsec  rise/fall  time.  As 
inferred  from  the  MOS  device  physics,  this  simulation  confirms 
the  expected  switching  times,  governed  exclusively  by  the 
effective  channel  length  (L^h)  > that  are  shortened  as  L^h 
decreases.  As  an  example,  for  y.j^  = 1.5  pim,  corresponding  to 
= 0.7  Jim,  the  circuit  is  expected  to  have  a turn-on  time 
of  13  nsec,  which  is  improved  approximately  28  % compared  with 
that  of  yjjj  = 4 ^m . Although,  this  kind  of  device  design 
predictability  is  crude,  it  reveals  its  potential  utility  in 
TCAD  applications. 

We  have  obtained  satisfying  simulation  results  in  case 
of  the  inductive-resistive  load  switching  circuit,  using  only 
one  lateral  DMOS  transistor.  We  now  present  a simulation  for 
a more  complex  circuit,  single-phase  full-bridge  inverter,  as 
shown  in  Figs.  4.12  and  4.13.  Typical  alternate  pulse  gener- 
ators are  used  as  gate  control  signal  to  obtain  AC  output.  The 
gate  pulse  amplitude  is  14  V with  0.5  msec  rise/fall  time  and 
the  input  supply  voltage  is  60  V.  In  the  circuit  diagram,  D^, 
D2,  D3,  and  D4  indicate  the  parasitic  source-to-drain  internal 
diodes  which  are  inherently  constructed  by  body-drift/ 
substrate-drift  p-n  junction  in  the  investigated  LDMOST 
structure  { see  Fig.  3.1  ).  These  diodes  across  the  power 
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Fig.  4.12  Circuit  diagram  of  single-phase  full-bridge  DC-to- 

AC  converter  and  input  file  for  LADISPICE  simula- 
tion. D^,  D2,  D3,  and  D4  indicate  the  source-to- 
drain  internal  parasitic  diodes  of  their  respec- 
tive LDMOSTs. 
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LDMOST  SINGLE-PHASE  FULL-BRIDGE  DC-TO-AC  CONVERTER 

Y1  3 4 11  NMODE  WZ=0.05  AE=0.25U  AC=0.55U  AEC=0.2U 

Y2  2 7 0 0 NMODE  WZ-0.05  AE=0.25U  AC^.SSU  AEC^.2U 

Y3  3 6 2 2 NMODE  \NZ=0.05  AE=0.25U  ACsO.SSU  AEC=:0.2U 

Y4  1 5 0 0 NMODE  WZ=0.05  AE=0.25U  AC=0.55U  AEC=0.2U 

RL2  815 
LL1  8 0.1M 
CL2  8 0.1M 

VDD  3 0 DC  60 

VG1  4 1 PULSE  ( 0 14  0 0.5M  0.5M  4M  10M ) 

VG2  7 0 PULSE  (0140  0.5M  0.5M  4M  1 0M ) 

VG3  6 2 PULSE  ( 0 14  5M  0.5M  0.5M  4M  10M  ) 

VG4  5 0 PULSE  ( 0 14  5M  0.5M  0.5M  4M  10M  ) 

.MODEL  NMODE  NLDMT 

+ (LCH=0.84E-4  LD1=1.0E-4  LD2=0.4E-4  LD3=2.6E-4  DELTAL=0.35E-4 
+ TOX=7.0E-6  TFOX=0.5E-4  YJB=1.7E-4  YJDP=3.2E-4  YEPI=7.2E-4 
+ NA0=3.0E17  NDD=1.5E16  NAS=5.0E15  NABEFF=8.0E17  VSAT=1.0E7 
+ THETA=1 .97E-1  UNCH0=4.0E2  UND=1.2E3  BETAN=1.5  BB=3.5 
+ TAUNE=0.15E-6  TAUHB=0.3E-6  TAUHC=0.3E-6  LAB=8.6E-4  LAC=23.6E-4 
+ VFBCH=-2.57  VFBDR=-0.2  PSIBD=0.7  FQD=0.8  KF0=0.1  KF1=3.0 
+ KF2»1.0  KF4=0.94  RBD=3.5  RDR=56.7  RSUB=62.1  ) 

.TRAN  0.1  M 14M  2M 

.PRINT  TRAN  V(4,1)  V(6,2)  V(1,2)  V(8,2) 

.PLOT  TRAN  V(4,1)  V(6,2)  V(1,2)  V(8,2) 

.OPTIONS  LIST  NODE  ACCT  LIMPTS=1001 

.WIDTH  IN=80  OUT=80 

.END 


Fig.  4.12 


continued 


Voltage  (V) 
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100.0  - 


50.0  F 


-100.0 


-50.0  - 


6.0  10.0 
Time  (msec) 


14.0 


Fig.  4.13  LADISPICE-simulated  output  voltage  waveforms  of 

single-phase  full-bridge  DC-to-AC  converter. 
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semiconductor  switch  are  frequently  used  to  clamp  the  turn- 
off inductive  load  current  [Mot89] . The  diodes  do  not  protect 
their  respective  switches  but  rather  the  complementary 
switches.  As  an  example,  in  the  configuration  of  Fig.  4.12, 
diodes  D^,  D2  protect  Y3,  Y4  and  D3,  D4  protect  Y^,  Y2.  To 
illustrate  this,  assume  Y^,  Y2  are  initially  conducting, 
causing  load  current  to  flow  from  Vqd  through  the  inductive 
load  to  the  ground.  When  Y^,  Y2  turn  off,  the  inductive  current 
will  continue  but  now  from  the  ground  through  D4  and  D3  to  the 
power  supply  Vdq.  Consequently,  the  fly-back  voltage  will  be 
clamped  to  In  Fig.  4.13,  the  simulation  starts  when  Y^, 
Y2  are  in  the  *on*  state  while  Y3,  Y4  are  in  the  *off*  state. 

( In  this  simulation  example,  there  has  been  a convergence 
difficulty  for  a high  value  of  load  elements.  ) This  simula- 
tion shows  the  capability  of  our  LDMOST  model  to  simulate  the 
parasitic  internal  diode.  It  is  more  pronounced  in  the 
reverse-recovery  transient  of  the  LDMOST  shown  in  Fig.  4.14. 
This  versatile  capability  is  important  because  the  circuit 
designer  can  use  these  internal  diodes  to  suppress  the  result- 
ing inductive  kick  from  exceeding  the  breakdown  voltage  of  the 
power  MOSFET. 

Finally,  the  job-time  performance  of  LADISPICE  is  illus- 
trated. Table  4.6  tabulates  the  total  number  of  iterations  and 
job- time,  performing  transient  simulations  on  a SUN  SPARC 
station-2 . We  believe  that  our  model  should  be  compared  with 


(VLU) 
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Fig.  4.14  LADISPICE-simulated  reverse-recovery  current 

transient  for  the  LDMOST.  The  supply  voltage, 
applied  through  a 20  £2  series  resistor,  changes 
from  -5  to  10  V with  a 0.6  |Xsec  rise  time  starting 
at  0.5  jlsec. 
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TABLE  4.6 

LADISPICE-1.1  JOB  TIME  PERFORMANCE 
{ on  SUN  SPARC  station-2  ) 


( job  times  in  sec  ) 


simulated  circuit 

iteration  no. 

job  time 

switching 

circuit 

switching 

1379 

30.95 

gate"^ 

611 

14.98 

turn-of 

899 

26.57 

turn-on'*"^'^ 

184 

7.08 

DC -AC  converter 

17911 

173.50 

+ with  Rq  = 500  ohms 
++  with  Ll  = 0.18  uH 
+++  with  Yjb  = 4.07  um 
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another  LDMOST  model  in  the  computation  time.  But,  for  lack 
of  an  appropriate  physical  or  subcircuit  model  in  conjuction 
with  the  investigated  LDMOST,  we  merely  list  the  job-time 
performance  of  our  semi -numerical  LADISPICE.  Overall, 
LADISPICE  tends  to  be  computationally  less  efficient.  This 
result  shows  intensive  conputation  time  when  the  semi -numer- 
ical method  is  used  for  device  modeling  instead  of 
conventional  analytical  method.  Moreover,  the  Runge-Kutta 
method  for  solving  a coiiplex  drift-region  model  equation 
requires  further  intensive  run-time.  However,  the  physical 

I 

and  semi -numerical  LADISPICE  provides  an  accurate  simulation 
in  circuit  performance  as  well  as  might  be  used  as  a reliable 
tool  for  advanced  power  IC  TCAD. 


4.7  Su 


The  lateral  DMOS  transistor  model  developed  in  previous 
chapters  has  been  successfully  irtplemented  into  a SPICE 
circuit  simulation  program,  generating  LADISPICE.  The  simula- 
tor SPICE2G.6  has  been  modified  and  a new  input  format  for  the 
LDMOST  has  been  implemented.  The  procedures  for  implementa- 
tion were  discussed  in  detail.  The  new  LDMOST  model  in  SPICE 
has  sufficient  flexibility  to  account  for  the  unique  charac- 
teristics of  the  LDMOST  operation,  while  still  retaining 
readily  extractable  model  parameters.  With  an  advantage  of 
directly  using  physical  and  process  parameters  in  its  device 
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model,  which  is  not  afforded  by  other  empirical  subcircuit 
SPICE  models,  the  predicted  trends  are  reasonable.  Further- 
more, LADISPICE  seems  quite  reliable.  In  our  representative 
simulation  exanples,  there  has  not  been  any  problem  { conver- 
gence, overflow,  underflow,  etc.  ),  even  without  using  the 
user  options  in  the  .OPTIONS  card  to  ease  convergence. 
Although  supporting  measured  data  are  not  available,  examples 
of  simulations  with  the  new  model  agree  well  qualitatively 
with  theoretical  prediction  and  other  experimental  investiga- 
tions. In  consequence,  LADISPICE  could  be  a powerful  tool  for 
circuit  simulation  in  power  (LDMOST)  IC  TCAD  applications.  In 
addition,  the  developed  LDMOST  model  is  portable  to  other 
SPICE-like  circuit  simulators . 


CHAPTER  5 

SUMMARY  AND  SUGGESTIONS  FOR  FUTURE  WORK 


In  this  dissertation,  a physically  based  predictive 
nonplanar-drift  lateral  DMOS  transistor  model  has  been 
presented  and  verified  with  experimental  measurements.  The 
charge-based  model  was  derived  based  on  a regional  approach 
involving  technological  model  parameters.  Different  from  the 
existing  subcircuit  power  device  models,  the  developed  model 
accounts  for  the  unique  device  structure  such  as  the  graded- 
channel  and  the  nonplanar-drift  region  as  well  as  the  unique 
device  characteristics  such  as  the  space-charge-limited 
current  flow  in  the  drift  region  and  the  inherent  parasitic 
BJT  action.  Even  though  a specific-oriented  DMOS  transistor 
has  been  investigated  and  modeled,  the  modeling  methodology 
and  the  developed  model  equations  in  this  work  may  be  appli- 
cable to  the  other  structures  of  the  DMOS  transistor  as  well, 
through  modification  of  the  models. 

The  developed  model  has  been  implemented  directly  in 
SPICE2G.6  source  code  to  create  LADISPICE.  The  final  version 
of  the  code,  LADISPICE-1 . 1,  evolved  from  this  dissertation, 
is  available  at  the  University  of  Florida.  With  an  advantage 
of  directly  using  physical/technological  parameters  in  its 
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device  model,  the  modified  SPICE  can  be  effectively  used  as  a 
design  tool  for  lateral  DMOS  transistor  IC's.  The  utility  and 
the  efficiency  of  LADISPICE-1 . 1 were  demonstrated  through 
representative  circuit  simulations.  The  developed  LDMOST 
model  is  also  portable  to  other  SPICE-like  circuit  simulators. 

Based  on  the  research  discussed  herein,  we  suggest  the 
following  topics  for  future  research  consideration. 

(1)  To  better  approximate  the  actual  2-D  doping  profile  in 
the  graded-channel,  the  functional  representation  should 
be  more  complex  than  a simple  1-D  exponential  approxima- 
tion. An  accurate  approximation  of  body-doping  density 
may  increase  the  accuracy  of  the  model  prediction, 
especially  at  low  V^s. 

(2)  Numerical  aspects  need  to  be  investigated  further  for 
potential  improvements  in  convergence  and  numerical 
efficiency.  As  an  example,  a system  of  differential 
equations  in  the  drift  region  can  be  replaced  with  an 
analytical  model  equation  with  a minor  damage  to  the 
physical  representation  of  device  operation.  It  may 
reduce  the  CPU-time  significantly  in  circuit  simulation. 

(3)  Modeling  of  carrier  impact  ionization  is  crucial  in 
characterization  of  breakdown  behavior.  Understanding 
the  relevant  physics  for  breakdown,  and  subsequently, 
the  incorporation  of  this  model  may  be  essential  to  the 
optim\im  design  of  the  LDMOST  device  and  circuit . 
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(4)  A drastic  study  of  tenperature  dependences  and  their 
incorporation  into  the  model  are  indispensable  for  power 
device  TCAD.  Most  of  device  parameters,  especially  the 
carrier  mobility,  will  be  affected,  and  the  device 
characteristics  will  be  significantly  changed  with 
regard  to  teirperature  variation. 

(5)  To  make  the  model  scalable,  3 -dimensional  geometry  and 
multi-dimensional  current  flows  need  to  be  investigated 
and  properly  accounted  for. 

(6)  For  the  model  to  be  more  applicable  in  all  regions  of 
lateral  DMOS  transistor  operations,  modeling  of  the 
subthreshold  region  is  recommended. 

(7)  For  a fully  automated  parameter  extraction,  a well 
defined  extraction  scheme  should  be  developed  and  incor- 
porated into  a parameter  extraction  tool  package. 

(8)  Based  on  the  model  discussed  in  this  dissertation,  the 
extension  of  the  model  to  p-channel  LDMOST  is  a worth- 


while task. 


APPENDIX  A 

LINKED-LIST  SPECIFICATIONS  OF  LATERAL  DMOS  TRANSISTOR 


A.l  Lateral  DMOS  Trai 


1 

o 

• • 

subckt  inform 

LOC  +00: 

next  pointer 

LOCV  +00: 

+01: 

LOCV  pointer 

+01: 

+02: 

nd 

+02: 

+03: 

ng 

+03: 

+04: 

ns 

+04: 

+05: 

nsub 

+05: 

+06: 

ne 

+06: 

+07: 

nb 

+07: 

+08: 

nc 

+08: 

+09: 

model  pointer 

+09: 

+10: 

off 

+10: 

+11: 

<unused> 

+12: 

<\inused> 

+13: 

(nd,ng) 

+14: 

(nd,ns) 

+15: 

(nd,nb) 

+16: 

(ng,nd) 

+17: 

(ng,ns) 

+18: 

(ns,nd) 

+19: 

(ns,ng) 

+20: 

(ns,ne) 

+21: 

(nsub,nc) 

+22: 

(ne,ns) 

+23: 

(ne,nb) 

+24: 

(ne,nc) 

+25: 

(nb,nd) 

+26: 

(nb,ne) 

+27: 

(nb,nc) 

+28: 

( nc , nsub ) 

+29: 

(nc,ne) 

+30: 

(nc,nb) 

+31: 

Lxi  offset 

+35: 

+32: 

(nd,nd) 

+36: 

+33: 

(ng,ng) 

+37: 

+34: 

(ns, ns) 

+38: 

sistor  (ID=15) 


device  name 

device  width 

IC-vds 

IC-vgs 

IC-vss 

<unused> 

<unused> 

<ianused> 

AE 

AC 

AEC 


(nsub,nsub) 

(ne,ne) 

(nb,nb) 

(nc,nc) 
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State  Variable  Table 


Lxi  +00 

vdso 

+35 

cgsogso 

+01 

vgso 

+36 

csidso 

+02 

<unused> 

+37 

csigso 

+03 

vebo 

+38 

cdidso 

+04 

vcbo 

+39 

cdigso 

+05 

cidso 

+40 

cbidso 

+06 

gdsdso 

+41 

cbigso 

+07 

gdsgso 

+42 

cgdodso 

+08 

cibfo 

+43 

cgdogso 

+09 

gbfebo 

+44 

cjeebo 

+10 

gbfcbo 

+45 

cjccbo 

+11 

cibro 

+46 

cdeebo 

+12 

gbrebo 

+47 

cdecbo 

+13 

gbrcbo 

+48 

cdcebo 

+14 

cicto 

+49 

cdccbo 

+15 

gctebo 

+50 

<unused> 

+16 

gctcbo 

+51 

<\inused> 

+17 

qgsolo 

+52 

<unused> 

+18 

cigsolo 

+19 

qsio 

+20 

cisio 

+21 

qdio 

+22 

cidio 

+23 

qbio 

+24 

cibio 

+25 

qgdolo 

+26 

cigdolo 

+27 

qjeo 

+28 

cijeo 

+29 

qjco 

+30 

ci  jco 

+31 

qdeo 

+32 

cideo 

+33 

qdco 

+34 

cidco 
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A. 2 Lateral  DMQS  Transistor  Model  (ID=25) 


-01;  subckt  inform 
LOG  +00:  next  pointer 
+01:  LOCV  pointer 
+02 ; model  type 


LOCV  +00: 

model  name 

+01: 

LCH 

+02: 

LDl 

+03: 

LD2 

+04: 

LD3 

+05: 

DELTAL 

+06: 

TOX 

+07: 

TFOX 

+08: 

YJB 

+09: 

YJDP 

+10: 

YEPI 

+11: 

NAO 

+12: 

NDD 

+13: 

NAS 

+14 : 

NABEFF 

+15: 

VSAT 

+16: 

THETA 

+17: 

UNCHO 

+18: 

UND 

+19: 

BETAN 

+20: 

BB 

+21: 

TAUNE 

+22: 

TAUHB 

+23; 

TAUHC 

+24: 

LAB 

+25: 

LAC 

+26: 

VFBCH 

+27: 

VFBDR 

+28: 

PSIBD 

+29: 

FQD 

+30: 

KFO 

+31: 

KFl 

+32: 

KF2 

+33: 

KF4 

+34: 

RBD 

+35: 

RDR 

+36: 

RSUB 

APPENDIX  B 

NODAL  ADMITTANCE  MATRIX  OF  LATERAL  DMOS  TRANSISTOR 


B.l  Pointer  Assignment  of  Nodal  Admittance  Matrix 


D 

G 

S 

SUB 

E 

B 

C 

D 

LOC+32 

LOC+13 

LOC+14 

LOC+15 

G 

LOC+16 

LOC+33 

LOC+17 

S 

LOC+18 

LOC+19 

LOC+34 

LOC+20 

SUB 

LOC+35 

LOC+21 

E 

LOC+22 

LOC+36 

LOC+23 

LOC+24 

B 

LOC+25 

LOC+26 

LOC+37 

LOC+27 

C 

LOC+28 

LOC+29 

LOC+30 

LOC+38 
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LOC+13 : 


LOC+14: 


LOC+15: 


LOC+16: 


LOC+17 : 


LOC+18: 


LOC+19 : 


LOC+20 : 


LOC+21: 


LOC+22 : 


LOC+23 : 


LOC+24: 


Nodal  Admittances  for  DC  and  Transient  Analysis 


^ where  i = GSOL,  SI,  DI,  BI,  GDOL,  JE,  JC,  DE,  DC  ) 


dl__ 

DS  DI  GDOL 


dl„- 

DS  DS  DX  DX  GDOL  GDOL 

5v  3v  5v  3v  5v 

^ DS  ^ GS  ^ DS  GS  ^ DS  GS 


1 


9'si 

*'bI,  ^'ODOL 

*''ds 

»'gsol  . 5'si  , , ®'di  . , 3'bi  , *'bi 

^^GDOL 

^^GDOL 

>''os 

'^''db  »''ob  »''ds  «''os  »''ds  S''os 

%s 

*''gs 

.^Isi.^'bi 
^DS  8Vos 

S'ds 

^^GSOLj^^SI 

S''gs 

^GS  »''gS  «''gS 

1 


BD 


R 


1 


SUB 


1 


BD 


dl„„ 

BF  BF  CT  CT  ^ JE  DE  DE 

dv„^  ” 3v ~ 3v ~ 3v dv ~ 3v "5v! 


EB 


CB 


EB 


CB 


EB 


EB 


CB 


^^BF  ^^CT  ^^DE 
^''CB  «''CB  *''CB 
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LOC+25: 


LOC+26: 


LOC+27 : 


LOC+28: 


LOC+29 : 


LOC+30 : 


LOC+32 : 


LOC+33 : 


LOC+34 : 


LOC+35: 


LOC+36: 


LOC+37 : 


LOC+38: 


1 


DR 


^^BF 

^'br  ^'je 

^^DE 

dl-,^ 

DC 

»''eb 

»''eb  3Veb 

^''eb 

*''eb 

^^BF 

«'br  , 

^^DE 

3^dc 

«''CB 

^CB  «''CB 

®''CB 

»''CB 

1 


SUB 


3^eb 


^^BR  ^^BR  JC 

dv„^  3v "^3v "^3v ^SvZZ  W. 


EB 


CB 


EB 


CB 


CB 


EB 


CB 


1 dl— — dX— T dl____ 

1 PS  DI  GDOL 


"dr  ^''ds  %s 


^^GSOL  ^^BI^^^GDOL 


av 


GS 


^^GS  ^^GS’^^GS 


1 ^ ^ ^ ^^GSOL  ^^BI 

R "'’dV ^dV ^dV dV dV dV dV. 


BD 


DS 


GS 


GS 


DS 


GS 


DS 


GS 


1 


SUB 


1 ^ ^^BF  ^ JE  ^ ^^DE 


^BD  ^^EB  ^^EB  ^^EB  ^^EB 


1 ^ ^^BF  ^ ^^BF  ^ ^^BR  ^ ^^BR  ^^JE  ^ ^^DE  ^ ^^DE  ^ ^^DC  ^ ^^DC 


*"dr  ^''eb  ■ ^""cb  ■ ^""eb  ■ ^""cb  ^""eb  3Vcb  ■ 3Veb  ■ 9Vcb  avgg  dV^g 


1 


/^BR  ^^CT  ^^JC,3^DC 
^suB  ^""cb’^^cb'^^cb  ^""cb 


APPENDIX  C 

NEW  MODEL  ROUTINES  FOR  LATERAL  DMOS  TRANSISTOR 


C.l  Subroutine  LDMOST.f 


SUBROUTINE  LEMOST 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

C 

C THIS  ROUTINE  PROCESSES  THE  LATERAL  DOUBLE-DIFFUSED  MOS  TRANSISTORS 

C FOR  DC  AND  TRANSIENT  ANALYSES. 

C 

C MODELED  BY  YEONBAE  CHUNG  7/15/93 

C PROGRAMMED  BY  YEONBAE  CHUNG  8/24/94 

C 

C SPICE  VERSION  2G.6  SCCSID=TABINF  3/15/83 

COMMON  /TABINF/  lELMNT, ISBCKT, NSBCKT, lUNSAT, NUNSAT, ITEMPS , NUMTEM, 

1 ISENS, NSENS, IFOUR, NFOUR, IFIELD, ICODE, IDELIM, ICOLUM, INSIZE, 

2 JUNODE , LSBKPT , NUMBKP , lORDER , JMNODE , lUR , lUC , ILC , ILR , NUMOFF , ISR , 

3 NMOFFC, ISEQ, ISEQl,NEQN,NODEVS,NDIAG,  ISWAP,  IEQUA,MACINS, LVNIMl, 

4 LXO , LVN, LYNL, LYU, LYL, LXl , LX2 , LX3 , LX4 , LX5 , LX6 , LX7 , LDO , LDl , LTD, 

5 IMYNL, IMVN, LCVN, NSNOD, NSMAT, NSVAL, ICNOD, ICMAT, ICVAL, 

6 LOUTPT, LPOL, LZER, IRSWPF, IRSWPR, ICSWPF, ICSWPR, IRPT, JCPT, 

7 I ROWNO , JCOLNO , NTTBR , NTTAR , LVNTMP 

C SPICE  VERSION  2G.6  SCCSID=CIRDAT  3/15/83 

COMMON  /CIRDAT/  LOCATE ( 50 ) ,JELCNT{ 50) ,NUNODS,NCNODS,NUMNOD,NSTOP, 

1 NUT,NLT,NXTRM,NDIST,NTLIN, IBR,NUMVS,NUMALT,NUMCYC 
C SPICE  VERSION  2G.6  SCCSID=STATUS  3/15/83 

COMMON  /STATUS/  OMEGA, TIME, DELTA,DELOLD(7) ,AG(7) ,VT,XNI,EGFET, 

1 XMU , SFACTR , MODE , MODEDC , ICALC , INITF , METHOD , lORD , MAXORD , NONCON , 

2 ITERNO, ITEMNO, NOSOLV,MODAC, IPIV, IVMFLG, IPOSTP, ISCRCH, lOFILE 
C SPICE  VERSION  2G.6  SCCSID=KNSTNT  3/15/83 

COMMON  /KNSTNT/  TWOPI,XLOG2,XLOG10,ROOT2, RAD, BOLTZ, CHARGE, CTOK, 

1 GMIN, RELTOL, ABSTOL, VNTOL, TRTOL, CHGTOL, EPSO , EPSSIL, EPSOX, 

2 PIVTOL, PIVREL 

C SPICE  VERSION  2G.6  SCCSID=FLAGS  3/15/83 

COMMON  /FLAGS/  IPRNTA, IPRNTL, IPRNTM, IPRNTN, IPRNTO, LIMTIM, LIMPTS, 

1 LVLCOD, LVLTIM, ITLl , ITL2 , ITL3 , ITL4 , ITL5 , ITL6 , IGOOF, NOGO, KEOF 
COMMON  /VNR/  VNNR, VDCNNR, VTRNNR, VACNNR, DELV, ODELV 
C SPICE  VERSION  2G.6  SCCSID=BLANK  3/15/83 
COMMON  /BLANK/  VALUE (200000) 

DOUBLE  PRECISION  NA0,NDD,NAS 
INTEGER  NODPLC(64) 

COMPLEX  CVALUE(32) 

EQU I VALENCE  ( VALUE ( 1 ) , NODPLC ( 1 ) , CVALUE ( 1 ) ) 
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c 

c 


c 

c 


DIMENSION  VDSO(l) ,VGSO(l) ,VEBO(l) ,VCBO(l) , 

1 CIDSO(l) ,GDSDSO(l) ,GDSGSO(l) ,CIBFO(l) ,GBFEBO(l) ,GBFCBO(l) 

2 CIBRO(l) ,GBREBO(l) ,GBRCBO(l) ,CICTO(l) ,GCTEBO(l) ,GCTCBO(l) 

3 QGSOLO(l) ,CIGSOLO(l) ,QSIO(l) ,CISIO(l)  ,QDIO(l)  ,CIDIO(l) , 

4 QBIO(l) ,CIBIO(l) ,QGDOLO(l) ,CIGDOLO(l) ,QJEO(l) ,CIJEO(l) , 

5 QJCO(l) ,CIJCO(l) ,QDEO(l) ,CIDEO(l) ,QDCO(l) ,CIDCO(l) , 

6 CGSOGSO ( 1 ) , CSIDSO { 1 ) , CSIGSO ( 1 ) , CDIDSO { 1 ) , CDIGSO ( 1 ) , 

7 CBIDSO ( 1 ) , CBIGSO { 1 ) , CGDODSO ( 1 ) , CGDOGSO ( 1 ) , CJEEBO ( 1 ) , 

8 C JCCBO ( 1 ) , CDEEBO ( 1 ) , CDECBO ( 1 ) , CDCEBO ( 1 ) , CDCCBO ( 1 ) 


EQUIVALENCE  ( VDSO ( 1 ) , VALUE ( 1 ) ) , 

1 (VEBO(l) ,VALUE(4) ) , 

EQUIVALENCE  (CIDSO(l) ,VALUE(6) ) , 

1 (GDSGSO(l) ,VALUE(8) ) , 

2 (GBFEBOd)  , VALUE (10) ) , 

3 (CIBRO(l) ,VALUE(12) ) , 

4 (GBRCBO(l) , VALUE (14) ) , 

5 (GCTEBO(l) , VALUE (16) ) , 
EQUIVALENCE  ( QGSOLO ( 1 ) , VALUE (18)), 

1 (QSIO(l) ,VALUE(20) ) , 

2 (QDIO(l) ,VALUE(22) ) , 

3 (QBIO(l) ,VALUE(24) ) , 

4 (QGDOLO(l) ,VALUE(26) ) , 

5 (QJEO(l) ,VALUE(28) ) , 

6 (QJCO(l) ,VALUE(30) ) , 

7 (QDEO(l) ,VALUE(32) ) , 

8 (QDCO(l) ,VALUE(34) ) , 
EQUIVALENCE  ( CGSOGSO ( 1 ) , VALUE (36)) 

1 (CSIGSO(l) ,VALUE(38) ) , 

2 (CDIGSO(l) ,VALUE(40) ) , 

3 (CBIGSO(l) ,VALUE(42) ) , 

4 (CGDOGSO ( 1) , VALUE (44) ) 

5 (CJCCBO(l) ,VALUE(46) ) , 

6 (CDECBO ( 1) ,VALUE(48) ) , 

7 (CDCCBO ( 1) ,VALUE( 50) ) 


(VGSO(l) ,VALUE(2) ) , 
(VCBO(l) ,VALUE(5) ) 
(GDSDSO(l) ,VALUE(7) ) , 
(CIBFO(l) ,VALUE(9) ) , 
(GBFCBO ( 1 ) , VALUE ( 11 ) ) , 
(GBREBO(l) .VALUE (13) ) , 
(CICTO(l) ,VALUE(15) ) , 
(GCTCBO(l) ,VALUE(17) ) 
(CIGSOLO ( 1 ) , VALUE ( 19 ) ) , 
(CISIO(l) ,VALUE(21) ) , 
(CIDIO(l) ,VALUE(23) ) , 
(CIBIO(l) ,VALUE(25) ) , 
(CIGDOLO(l) ,VALUE(27) ) , 
(CIJEO(l) ,VALUE(29) ) , 
(CIJCO(l) ,VALUE(31) ) , 
(CIDEO(l) ,VALUE(33) ) , 
(CIDCO(l) ,VALUE(35) ) 
(CSIDSO(l) ,VALUE(37) ) , 
(CDIDSO(l) ,VALUE(39) ) , 
(CBIDSO(l) ,VALUE(41) ) , 
(CGDODSO ( 1 ) , VALUE ( 43 ) ) , 
(CJEEBO ( 1 ) , VALUE ( 45 ) ) , 
(CDEEBO ( 1 ) , VALUE ( 47 ) ) , 
(CDCEBO(l) ,VALUE(49) ) , 


/ 

$ 


LOC=LOCATE ( 15 ) 

10  IF  ( (LOC.EQ.O) .OR. (NODPLC(LOC+39) .NE.O) ) RETURN 
LOCV=NODPLC ( LOC+ 1 ) 

ND=NODPLC (LOC+2 ) 

NG=NODPLC (LOC+3 ) 

NS=NODPLC (LOC+4 ) 

NSUB=NODPLC (LOC+5 ) 

NE=NODPLC ( LOC+  6 ) 

NB=NODPLC (LOC+7 ) 

NC=NODPLC (LOC+8 ) 

LOCM=NODPLC (LOC+9 ) 

IOFF=NODPLC ( LOC+ 1 0 ) 

TYPE=NODPLC (LOCM+2 ) 

LOCM=NODPLC ( LOCM+ 1 ) 

LOCT=NODPLC ( LOC+3 1 ) 

LOCTO=LXO+LOCT 
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L0CT1=LX1+L0CT 

LOCT2=LX2+LOCT 

C 

C MODEL  PARAMETERS 
C 

WZ  = VALUE (LOCV+1) 

AE  = VALUE (LOCV+8) 

AC  = VALUE (LOCV+9) 

DELTAL  = VALUE(LOCM+5) 

TOX  = VALUE (LOCM+6) 

NAO  = VALUE (LOCM+ 11) 

NDD  = VALUE (LOCM+ 12) 

NAS  = VALUE (LOCM+ 13) 

VFBCH  = VALUE(LOCM+26) 

PSIBD  = VALUE(LOCM+28) 

GBD  = VALUE(LOCM+34) 

GDR  = VALUE(LOCM+35) 

GSUB  = VALUE(LOCM+36) 

COX  = EPSOX/TOX 

PIBO  = VT*DLOG(NAO/XNI) 

VTH  = VFBCH  + 2.0DO*PIBO 

& + DSQRT(4.0D0*CHARGE*EPSSIL*NA0*PIB0)/COX 

PSISUB  = VT*DLOG(NDD*NAS/(XNI*XNI) ) 

VDINIT  = l.ODO 
VGINIT  = VTH  + l.ODO 
VCRIT  = 0.6D0 
C 

C INITIALIZATION 
C 

ICHECK=1 

GOTO  (100,20,30,50,60,70),  INITF 
C . . INITIALIZE  JUNCTION  VOLTAGES 
20  IF  (lOFF.NE.O)  GO  TO  40 
VDS=TYPE*VALUE (LOCV+2 ) 

VGS=TYPE*VALUE ( LOCV+3 ) 

VEB=-VDS 

VCB=TYPE*VALUE (LOCV+4 ) -VDS 
IF  (VDS.NE.O.ODO)  GO  TO  300 

IF  (VGS.NE.O.ODO)  GO  TO  300 

IF  (VEB.NE.O.ODO)  GO  TO  300 

IF  (VCB.NE.O.ODO)  GO  TO  300 

IF  ( (MODE.EQ.l) .AND. (MODEDC.EQ.2) .AND. (NOSOLV.NE.O) ) GO  TO  300 

VDS=VDINIT 

VGS=VGINIT 

VEB=-VDS 

VCB=-VDS 

GO  TO  300 

C . . CONVERGE  WITH  ' OFF ' DEVICES  HELD  OFF 
30  IF  (lOFF.EQ.O)  GO  TO  100 
40  VDS=0.0D0 
VGS=0.0D0 
VEB=0 . ODO 
VCB=0 . ODO 
GO  TO  300 
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C.. STORE  SMALL-SIGNAL  PARAMETERS 
50  VDS=VDSO(LOCTO) 

VGS=VGSO(LOCTO) 

VEB=VEBO(LOCTO) 

VCB=VCBO(LOCTO) 

GO  TO  300 

C.. FIRST  POINT  IN  TRANSIENT  ANALYSIS 
60  VDS=VDSO ( LOCTl ) 

VGS=VGSO(LOCTl) 

VEB=VEBO( LOCTl) 

VCB=VCBO( LOCTl) 

GO  TO  300 

C,, FIRST  ITERATION  OF  A SWEEP  POINT.  PREDICT  JUNCTION  VOLTAGES 
7 0 XFACT=DELTA/ DELOLD ( 2 ) 

VDSO ( LOCTO ) =VDSO ( LOCTl ) 

VDS= ( 1 . ODO+XFACT) *VDSO (LOCTl ) -XFACT*VDSO (LOCT2 ) 

VGSO ( LOCTO ) =VGSO { LOCTl ) 

VGS= ( 1 . ODO+XFACT) *VGSO (LOCTl ) -XFACT*VGSO (LOCT2 ) 

VEBO ( LOCTO ) =VEBO ( LOCTl ) 

VEB= ( 1 . ODO+XFACT) *VEBO ( LOCTl ) -XFACT*VEBO ( LOCT2 ) 

VCBO ( LOCTO ) =VCBO ( LOCTl ) 

VCB= ( 1 . ODO+XFACT) *VCBO (LOCTl ) -XFACT*VCBO (LOCT2 ) 

CIDSO (LOCTO ) =CIDSO (LOCTl) 

GDSDSO ( LOCTO ) =GDSDSO ( LOCTl ) 

GDSGSO ( LOCTO ) =GDSGSO ( LOCTl ) 

CIBFO ( LOCTO ) =C IBFO ( LOCTl ) 

GBFEBO ( LOCTO ) =GBFEBO ( LOCTl ) 

GBFCBO ( LOCTO ) =GBFCBO ( LOCTl ) 

C IBRO ( LOCTO ) =C IBRO ( LOCTl ) 

GBREBO ( LOCTO ) =GBREBO ( LOCTl ) 

GBRCBO ( LOCTO ) =GBRCBO ( LOCTl ) 

C ICTO ( LOCTO ) =C ICTO ( LOCTl ) 

GCTEBO ( LOCTO ) =GCTEBO ( LOCTl ) 

GCTCBO ( LOCTO ) =GCTCBO ( LOCTl ) 

CGSOGSO ( LOCTO ) =CGSOGSO ( LOCTl ) 

CSIDSO (LOCTO ) =CSIDSO (LOCTl ) 

CS IGSO ( LOCTO ) =CS IGSO ( LOCTl ) 

CDIDSO ( LOCTO ) =CDIDSO ( LOCTl ) 

CDIGSO ( LOCTO ) =CDIGSO ( LOCTl ) 

CBIDSO (LOCTO ) =CBIDSO (LOCTl ) 

CBIGSO (LOCTO ) =CBIGSO (LOCTl ) 

CGDODSO ( LOCTO ) =CGDODSO ( LOCTl ) 

CGDOGSO ( LOCTO ) =CGDOGSO ( LOCTl ) 

C JEEBO ( LOCTO ) =C JEEBO ( LOCTl ) 

C JCCBO ( LOCTO ) =C JCCBO ( LOCTl ) 

CDEEBO ( LOCTO ) =CDEEBO ( LOCTl ) 

CDECBO ( LOCTO ) =CDECBO ( LOCTl ) 

CDCEBO ( LOCTO ) =CDCEBO ( LOCTl ) 

CDCCBO ( LOCTO ) =CDCCBO ( LOCTl ) 

GO  TO  110 

C 

C COMPUTE  NEW  NONLINEAR  BRANCH  VOLTAGES 

C 

100  VDS=TYPE* (VALUE (LVNIMl+ND) -VALUE (LVNIMl+NS) ) 
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VGS=TYPE* (VALUE (LVNIMl+NG) -VALUE (LVNIMl+NS) ) 

VEB=TYPE* (VALUE (LVNIMl+NE) -VALUE (LVNIMl+NB) ) 

VCB=TYPE* (VALUE (LVNIMl+NC) -VALUE (LVNIMl+NB) ) 

110  DELVDS=VDS-VDSO(LOCTO) 

DELVGS=VGS-VGSO (LOCTO ) 

DELVEB=VEB-VEBO ( LOCTO ) 

DELVCB=VCB-VCBO ( LOCTO ) 

CIDSHAT=CIDSO(LOCTO) +GDSDSO(LOCTO) *DELVDS+GDSGSO(LOCTO) *DELVGS 
CIBFHAT=CIBFO (LOCTO ) +GBFEBO (LOCTO ) *DELVEB+GBFCBO ( LOCTO ) *DELVCB 
CIBRHAT=CIBRO (LOCTO ) +GBREBO (LOCTO ) *DELVEB+GBRCBO (LOCTO ) *DELVCB 
CICTHAT=CICTO ( LOCTO ) +GCTEBO ( LOCTO ) *DELVEB+GCTCBO ( LOCTO ) *DELVCB 
C 

C BYPASS  IF  SOLUTION  HAS  NOT  CHANGED 
C 

IF  (INITF.EQ.6)  GO  TO  200 

TOL=RELTOL*DMAXl(DABS(VDS) , DABS (VDSO (LOCTO) ) ) +VNTOL 
IF  (DABS(DELVDS) .GE.TOL)  GO  TO  200 

T0L=RELT0L*I»1AX1 (DABS (VGS) , DABS (VGSO (LOCTO ) ) ) +VNTOL 
IF  (DABS (DELVGS) .GE.TOL)  GO  TO  200 

TOL=RELTOL*E»lAXl (DABS (VEB) , DABS (VEBO (LOCTO ) ) ) +VNTOL 
IF  (DABS (DELVEB) .GE.TOL)  GO  TO  200 

TOL=RELTOL*EMAXl (DABS (VCB) , DABS (VCBO (LOCTO ) ) ) +VNTOL 
IF  (DABS (DELVCB) .GE.TOL)  GO  TO  200 

TOL=RELTOL*I»lAXl (DABS (CIDSHAT) , DABS (CIDSO (LOCTO ) ) ) +ABSTOL 
IF  (DABS(CIDSHAT-CIDSO(LOCTO) ) .GE.TOL)  GO  TO  200 
TOL=RELTOL*I»lAXl (DABS (CIBFHAT) , DABS (CIBFO (LOCTO ) ) ) +ABSTOL 
IF  (DABS (CIBFHAT-CIBFO (LOCTO )) .GE.TOL)  GO  TO  200 
TOL=RELTOL*E»lAXl (DABS (CIBRHAT) , DABS (CIBRO (LOCTO ) ) ) +ABSTOL 
IF  (DABS (CIBRHAT-CIBRO (LOCTO )) .GE.TOL)  GO  TO  200 
TOL=RELTOL*I»lAXl (DABS (CICTHAT) , DABS (CICTO (LOCTO ) ) ) +ABSTOL 
IF  ( DABS (CICTHAT-CICTO( LOCTO )) .GE.TOL)  GO  TO  200 
VDS=VDSO( LOCTO) 

VGS =VGSO( LOCTO) 

VEB=VEBO ( LOCTO ) 

VCB=VCBO( LOCTO) 

CIDS=CIDSO (LOCTO ) 

GDSDS=GDSDSO ( LOCTO ) 

GDSGS=GDSGSO ( LOCTO ) 

CIBF=CIBFO (LOCTO ) 

GBFEB=GBFEBO ( LOCTO ) 

GBFCB=GBFCBO ( LOCTO ) 

CIBR=CIBRO (LOCTO ) 

GBREB=GBREBO (LOCTO ) 

GBRCB=GBRCBO ( LOCTO ) 

CICT=CICTO (LOCTO ) 

GCTEB=GCTEBO (LOCTO ) 

GCTCB=GCTCBO ( LOCTO ) 

IF  (MODE.NE.l)  GO  TO  140 

IF  ( (MODEDC.EQ.2) .AND. (NOSOLV.NE.O) ) GO  TO  140 
GO  TO  850 
C 

C CONVERT  CAPACITANCES  TO  EQUIVALENT  CONDUCTANCES 
C 

140  GCGSOGS=AG ( 1 ) *CGSOGSO ( LOCTO ) 
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GCSIDS=AG(1) *CSIDSO(LOCTO) 

GCSIGS=AG ( 1 ) *CSIGSO (LOCTO ) 

GCDIDS=AG ( 1 ) *CDIDSO ( LOCTO ) 

GCDIGS=AG ( 1 ) *CDIGSO ( LOCTO ) 

GCBIDS=AG ( 1 ) *CBIDSO (LOCTO ) 

GCBIGS=AG(1) *CBIGSO (LOCTO) 

GCGDODS=AG ( 1 ) *CGDODSO ( LOCTO ) 

GCGDOGS=AG ( 1 ) *CGDOGSO ( LOCTO ) 

GC JEEB=AG ( 1 ) *C JEEBO ( LOCTO ) 

GC JCCB=AG ( 1 ) *C JCCBO ( LOCTO ) 

GCDEEB=AG ( 1 ) *CDEEBO ( LOCTO ) 

GCDECB=AG ( 1 ) *CDECBO (LOCTO ) 

GCDCEB=AG ( 1 ) *CDCEBO ( LOCTO ) 

GCDCCB=AG ( 1 ) *CDCCBO ( LOCTO ) 

GO  TO  860 
C 

C LIMIT  NONLINEAR  BRANCH  VOLTAGES 
C 

200  ICHK1=1 

CALL  FETLIM ( VGS , VGSO ( LOCTO ) , VTH ) 

IF  ( VDSO ( LOCTO ) .LT.O.ODO)  GO  TO  205 
CALL  LIMVDS ( VDS , VDSO ( LOCTO ) ) 

GO  TO  210 

205  CALL  LIMVDS ( -VDS , -VDSO ( LOCTO ) ) 

210  CALL  PNJLIM (VEB,VEBO (LOCTO ) ,VT,VCRIT, ICHECK) 

CALL  PNJLIM (VCB, VCBO (LOCTO ) , VT, VCRIT, ICHKl ) 

IF  ( ICHKl, EQ.l)  ICHECK=1 
C 

C DETERMINE  CURRENTS  AND  CONDUCTANCES  FOR  STATIC  CONDITIONS 
C 

300  CALL  LDTMOS(VDS,VGS,LOC,CIDS,QSI,QDI,QBI,QGDOL) 

CALL  LDTMOS ( VDS+DELV, VGS , LOC, CIDSVD, QSIVD, QDIVD, QBIVD, QGDOLVD) 
CALL  LDTMOS (VDS, VGS+DELV, LOC, CIDSVG, QSIVG, QDIVG, QBIVG, QGDOLVG) 
CALL  LDTBJT(VEB, VCB, LOC, CIBF, CIBR, CICT, QDE, QDC) 

CALL  LDTBJT ( VEB+DELV, VCB , LOC , CIBFVE , CIBRVE , CICTVE , QDEVE , QDCVE ) 
CALL  LDTBJT (VEB,VCB+DELV, LOC, CIBFVC, CIBRVC, CICTVC, QDEVC, QDCVC) 
C . . CONDUCTANCES  BY  NUMERICAL  DIFFERENTIATION 
C (FORWARD-DIFFERENCE  FORMULA) 

GDSDS= (CIDSVD-CIDS) *ODELV 
GDSGS= (CIDSVG-CIDS) *ODELV 
GBFEB= (CIBFVE-CIBF) *ODELV 
GBFCB= (CIBFVC-CIBF) *ODELV 
GBREB= (CIBRVE-CIBR) *ODELV 
GBRCB= (CIBRVC-CIBR) *ODELV 
GCTEB= (CICTVE-CICT) *ODELV 
GCTCB= (CICTVC-CICT) *ODELV 
C 

C BRANCH  ON  ANALYSIS  MODE 
C 

IF  (MODE.NE.l)  GO  TO  500 

IF  ( (MODEDC.EQ.2) .AND. (NOSOLV.NE.O) ) GO  TO  500 
IF  (INITF.EQ.4)  GO  TO  500 
GO  TO  700 
C 
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C DETERMINE  CHARGES  AND  CAPACITANCES  FOR  STATIC  CONDITIONS 
C 

500  QSIO(LOCTO)=QSI 
QDIO(LOCTO)=QDI 
QBIO(LOCTO)=QBI 
QGDOLO (LOCTO ) =QGDOL 
QDEO(LOCTO)=QDE 
QDCO( LOCTO )=QDC 

C. .CAPACITANCES  BY  NUMERICAL  DIFFERENTIATION 
C (FORWARD-DIFFERENCE  FORMULA) 

CSIDS= (QSIVD-QSI) *ODELV 
CSIGS= (QSIVG-QSI) *ODELV 
CDIDS= (QDIVD-QDI) *ODELV 
CDIGS= (QDIVG-QDI) *ODELV 
CBIDS= (QBIVD-QBI) *ODELV 
CBIGS= (QBIVG-QBI) *ODELV 
CGDODS= (QGDOLVD-QGDOL) *ODELV 
CGDOGS= (QGDOLVG-QGDOL) *ODELV 
CDEEB= (QDEVE-QDE) *ODELV 
CDECB= (QDEVC-QDE) *ODELV 
CDCEBs (QDCVE-QDC) *ODELV 
CDCCB= (QDCVC-QDC) *ODELV 
C 

QGSOLO (LOCTO ) =WZ*DELTAL*COX*VGS 
CGSOGS=WZ*DELTAL*COX 
IF  (VEB.LT.PSIBD)  THEN 

QJE1=DSQRT(2 . ODO*CHARGE*EPSSIL*PSIBD*NDD) 

Q JE2  =DSQRT ( 1 . ODO -VEB/ PS IBD ) 

QJEO (LOCTO ) =AE*QJE1*QJE2 
CJEEB=-AE*QJE1/ (2 . 0D0*PSIBD*QJE2 ) 

ELSE 

QJEO (LOCTO ) =0 . ODO 
CJEEB=0,0D0 
ENDIF 

IF  (VCB.LT.PSISUB)  THEN 

QJC1=DSQRT(2.0D0*CHARGE*EPSSIL*PSISUB*NDD*NAS/ (NDD+NAS) ) 
QJC2=DSQRT( 1 . ODO-VCB/PSISUB) 

QJCO (LOCTO ) =AC*QJC1*QJC2 
CJCCB=-AC*QJC1/ (2.0D0*PSISUB*QJC2) 

ELSE 

QJCO (LOCTO) =0.0 DO 
CJCCB=0.0D0 
ENDIF 
C 

C STORE  SMALL -SIGNAL  PARAMETERS 
C 

580  IF  ( (MODE.EQ.l) .AND. (MODEDC.EQ.2) .AND. (NOSOLV.NE.O) ) GO  TO  700 
IF  (INITF.NE.4)  GO  TO  600 
VALUE ( LOCTO +35) =CGSOGS 
VALUE ( LOCTO +36)=CSIDS 
VALUE (LOCTO+37 ) =CSIGS 
VALUE ( LOCTO+3 8 ) =CDIDS 
VALUE (LOCTO+39 ) =CDIGS 
VALUE (LOCTO+40 ) =CBIDS 
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VALUE (LOCTO+41 ) =CBIGS 
VALUE (LOCTO+42 ) =CGDODS 
VALUE (LOCTO+ 43 ) =CGDOGS 
VALUE ( LOCTO +44) =C JEEB 
VALUE ( LOCTO + 4 5 ) =C JCCB 
VALUE (LOCTO+46) =CDEEB 
VALUE (LOCTO+47 ) =CDECB 
VALUE ( LOCTO + 4 8 ) =CDCEB 
VALUE (LOCTO+49 ) =CDCCB 
GO  TO  1000 
C 

C TRANSIENT  ANALYSIS 
C 


600  IF  (INITF.NE.5)  GO  TO  610 
QGSOLO ( LOCTl ) =QGSOLO ( LOCTO ) 

QS 10 ( LOCTl ) =QS 10 ( LOCTO ) 

QDIO ( LOCTl ) =QDIO ( LOCTO ) 

QBIO { LOCTl ) =QBIO ( LOCTO ) 

QGDOLO ( LOCTl ) =QGDOLO ( LOCTO ) 

QJEO { LOCTl ) =Q JEO ( LOCTO ) 

QJCO ( LOCTl ) =Q JCO ( LOCTO ) 

QDEO { LOCTl ) =QDEO ( LOCTO ) 

QDCO ( LOCTl ) =QDCO ( LOCTO ) 

C . . INTEGRATE  CHARGES 

610  CALL  INTGR8(GEQ,CEQ,O.ODO,LOCT+17) 
CALL  INTGR8 ( GEQ , CEQ , 0 . ODO , LOCT+ 1 9 ) 
CALL  INTGR8(GEQ,CEQ,O.ODO,LOCT+21) 
CALL  INTGR8(GEQ,CEQ,O.ODO,LOCT+23) 
CALL  INTGR8(GEQ,CEQ,O.ODO,LOCT+25) 
CALL  INTGR8(GEQ,CEQ,O.ODO,LOCT+27) 
CALL  INTGR8(GEQ,CEQ,O.ODO,LOCT+29) 
CALL  INTGR8(GEQ,CEQ,O.ODO,LOCT+31) 
CALL  INTGR8(GEQ,CEQ,O.ODO,LOCT+33) 
C . , CALCULATE  EQUIVALENT  CONDUCTANCES 
GCGSOGS=AG ( 1 ) *CGSOGS 
GCSIDS=AG ( 1 ) *CSIDS 
GCSIGS=AG ( 1 ) *CSIGS 
GCDIDS=AG ( 1 ) *CDIDS 
GCDIGS=AG ( 1 ) *CDIGS 
GCBIDS=AG ( 1 ) *CBIDS 
GCBIGS=AG ( 1 ) *CBIGS 
GCGDODS=AG ( 1 ) *CGDODS 
GCGDOGS=AG ( 1 ) *CGDOGS 
GC JEEB=AG ( 1 ) *C JEEB 
GC JCCB=AG ( 1 ) *C JCCB 
GCDEEB=AG ( 1 ) *CDEEB 
GCDECB=AG ( 1 ) *CDECB 
GCDCEB=AG ( 1 ) *CDCEB 
GCDCCB=AG ( 1 ) *CDCCB 
C 


IF  (INITF.EQ.5)  THEN 

CIGSOLO (LOCTl ) =CIGSOLO (LOCTO ) 
C IS 10 ( LOCTl )=C IS 10 (LOCTO) 

C IDIO ( LOCTl ) =C IDIO ( LOCTO ) 
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CIBIO (LCXTTl) =CIBIO (LOCTO ) 

CIGDOLO (LOCTl ) =CIGDOLO (LOCTO ) 

CIJEO (LOCTl ) =CIJEO (LOCTO ) 

CIJCO (LOCTl ) =CIJCO (LOCTO ) 

CIDEO (LOCTl ) =CIDEO (LOCTO ) 

CIDCO (LOCTl ) =CIDCO (LOCTO ) 

ENDIF 

C 

C CHECK  CONVERGENCE 
C 

700  IF  (INITF.NE.3)  GO  TO  710 
IF  (lOFF.EQ.O)  GO  TO  710 
GO  TO  750 

710  IF  ( ICHECK . EQ . 1 ) GO  TO  720 

TOL=RELTOL*I»lAXl (DABS (CIDSHAT) , DABS (CIDS) ) +ABSTOL 
IF  (DABS(CIDSHAT-CIDS) .GE.TOL)  GO  TO  720 
TOL=RELTOL*I»lAXl (DABS (CIBFHAT) , DABS (CIBF) ) +ABSTOL 
IF  (DABS (CIBFHAT-CIBF) .GE.TOL)  GO  TO  720 
TOL=RELTOL*E»lAXl (DABS (CIBRHAT) , DABS (CIBR) ) +ABSTOL 
IF  (DABS (CIBRHAT-CIBR) .GE.TOL)  GO  TO  720 
TOL=RELTOL*EMAXl (DABS (CICTHAT) , DABS (CICT) ) +ABSTOL 
IF  (DABS (CICTHAT-CICT) .GE.TOL)  GO  TO  720 

720  NONCON=NONCON+ 1 

750  VDSO( LOCTO )=VDS 
VGSO( LOCTO )=VGS 
VEBO ( LOCTO ) = VEB 
VCBO( LOCTO )=VCB 
CIDSO (LOCTO ) =CIDS 
GDSDSO ( LOCTO ) =GDSDS 
GDSGSO (LOCTO ) =GDSGS 
CIBFO (LOCTO ) =CIBF 
GBFEBO ( LOCTO ) =GBFEB 
GBFCBO (LOCTO ) =GBFCB 
CIBRO (LOCTO ) =CIBR 
GBREBO (LOCTO ) =GBREB 
GBRCBO (LOCTO ) =GBRCB 
CICTO (LOCTO ) =CICT 
GCTEBO ( LOCTO ) =GCTEB 
GCTCBO (LOCTO ) =GCTCB 
IF  (MODE.NE.l)  GO  TO  800 

IF  ( (MODEDC.EQ.2) .AND. (NOSOLV.NE.O) ) GO  TO  800 
GO  TO  850 

800  CGSOGSO ( LOCTO ) =CGSOGS 
CSIDSO (LOCTO ) =CS IDS 
CSIGSO (LOCTO ) =CSIGS 
CDIDSO (LOCTO ) =CDIDS 
CDIGSO (LOCTO ) =CDIGS 
CBIDSO (LOCTO ) =CBIDS 
CBIGSO (LOCTO ) =CBIGS 
CGDODSO ( LOCTO ) =CGDODS 
CGDOGSO (LOCTO ) =CGDOGS 
C JEEBO ( LOCTO ) =C JEEB 
CJCCBO (LOCTO ) =CJCCB 
CDEEBO (LOCTO ) =CDEEB 


o o 
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CDECBO ( LOCTO ) =CDECB 
CDCEBO ( LOCTO ) =CDCEB 
CDCCBO ( LOCTO ) =CDCCB 
GO  TO  860 
C 

C SET  EQUIVALENT  CHARGE  CURRENTS  AND  CONDUCTANCES  TO  ZERO 
C FOR  DC  ANALYSIS 
C 

850  CIGSOLEQ=0 . ODO 
CISIEQ=0.0D0 
CIDIEQ=0.0D0 
CIBIEQ=0.0D0 
CIGDOLEQ=0 . ODO 
CIJEEQ=0.0D0 
CIJCEQ=0.0D0 
CIDEEQ=0.0D0 
CIDCEQ=0 . ODO 
GCGSOGS=0 . ODO 
GCSIDS=0.0D0 
GCSIGS=0.0D0 
GCDIDS=0 . ODO 
GCDIGS=0 . ODO 
GCBIDS=0.0D0 
GCBIGS=0.0D0 
GCGDODS=0 . ODO 
GCGDOGS=0 . ODO 
GCJEEB=0 . ODO 
GCJCCB=0 . ODO 
GCDEEB=0 . ODO 
GCDECB=0 . ODO 
GCDCEB=0 . ODO 
GCDCCB=0 . ODO 
GO  TO  900 
C 

C EVALUATE  EQUIVALENT  CHARGE  CURRENTS 
C 

860  C IGSOLEQ=C IGSOLO ( LOCTO ) -GCGSOGS  * VGS 

C I S IEQ=C I S 10 ( LOCTO ) -GCS IDS * VDS -GCS IGS * VGS 
CIDIEQ=CIDIO (LOCTO ) -GCDIDS*VDS-GCDIGS*VGS 
CIBIEQ=CIBIO (LOCTO )-GCBIDS*VDS-GCBIGS*VGS 
CIGDOLEQ=CIGDOLO (LOCTO ) -GCGDODS*VDS-GCGDOGS*VGS 
CIJEEQ=CIJEO (LOCTO) -GCJEEB*VEB 
CIJCEQ=CIJCO (LOCTO ) -GCJCCB*VCB 
CIDEEQ=CIDEO (LOCTO) -GCDEEB*VEB-GCDECB*VCB 
C IDCEQ=C IDCO ( LOCTO ) -GCDCEB * VEB -GCDCCB  * VCB 

EVALUATE  EQUIVALENT  DC  CURRENTS 
C 

900  C IDSEQ=C IDS -GDSDS  *VDS -GDSGS  *VGS 
C IBFEQ=C  IBF-GBFEB*VEB-GBFCB  *VCB 
CIBREQ=CIBR-GBREB*VEB-GBRCB*VCB 
C ICTEQ=C ICT-GCTEB  * VEB -GCTCB  * VCB 
C 

C LOAD  CURRENT  EXCITATION  VECTOR 
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c 


c 

c 

c 


VALUE (LVN+ND) 
VALUE (LVN+NG) 

& 

VALUE (LVN+NS) 
VALUE (LVN+NE) 
VALUE (LVN+NB) 

& 

VALUE (LVN+NC) 


=VALUE (LVN+ND) -CIDSEQ-CIDIEQ+CIGDOLEQ 
=VALUE (LVN+NG) -CIGSOLEQ+CISIEQ+CIDIEQ+CIBIEQ 

-CIGDOLEQ 

=VALUE (LVN+NS) +CIDSEQ+CIGSOLEQ-CISIEQ-CIBIEQ 
=VALUE (LVN+NE) -CIBFEQ-CICTEQ+CIJEEQ-CIDEEQ 
=VALUE (LVN+NB) +CIBFEQ+CIBREQ-CIJEEQ-CIJCEQ 

+CIDEEQ+CIDCEQ 

=VALUE (LVN+NC) -CIBREQ+CICTEQ+CIJCEQ-CIDCEQ 


LOAD  Y MATRIX 


LOCY=LVN+NODPLC (LOC+32 ) 
VALUE ( LOCY ) =VALUE ( LOCY ) 
LOCY=LVN+NODPLC ( LOC+3  3 ) 
VALUE ( LOCY ) = VALUE ( LOCY ) 
LOCY=LVN+NODPLC ( LOC+3  4 ) 
VALUE ( LOCY ) = VALUE ( LOCY ) 
Sc 

LOCY=LVN+NODPLC ( LOC+3  5 ) 
VALUE ( LOCY ) =VALUE ( LOCY ) 
LOCY=LVN+NODPLC (LOC+3  6 ) 
VALUE ( LOCY ) =VALUE ( LOCY ) 
LOCY=LVN+NODPLC (LOC+37 ) 
VALUE ( LOCY ) =VALUE ( LOCY ) 
Sc 

LOCY=LVN+NODPLC  ( LOC+3  8 ) 
VALUE { LOCY ) =VALUE { LOCY ) 
LOCY=LVN+NODPLC (LOC+13 ) 
VALUE ( LOCY ) =VALUE ( LOCY ) 
LOCY=LVN+NODPLC ( LOC+ 1 4 ) 
VALUE ( LOCY ) =VALUE ( LOCY ) 
Sc 

LOCY=LVN+NODPLC  ( LOC+ 1 5 ) 
VALUE ( LOCY ) =VALUE { LOCY ) 
LOCY=LVN+NODPLC ( LOC+ 1 6 ) 
VALUE ( LOCY ) =VALUE ( LOCY ) 
LOCY=LVN+NODPLC (LOC+17 ) 
VALUE ( LOCY ) =VALUE ( LOCY ) 
Sc 

LOCY=LVN+NODPLC  (LOC+18 ) 
VALUE ( LOCY ) =VALUE ( LOCY ) 
LOCY=LVN+NODPLC (LOC+19 ) 
VALUE ( LOCY ) =VALUE ( LOCY ) 
LOCY=LVN+NODPLC (LOC+20 ) 
VALUE { LOCY ) =VALUE { LOCY ) 
LOCY=LVN+NODPLC ( LOC+2 1 ) 
VALUE { LOCY ) =VALUE { LOCY ) 
LOCY=LVN+NODPLC {LOC+22 ) 
VALUE ( LOCY ) =VALUE { LOCY ) 
LOCY=LVN+NODPLC {LOC+23 ) 
VALUE ( LOCY ) =VALUE ( LOCY ) 
Sc 

LOCY=LVN+NODPLC  (LOC+24 ) 


+GDR+GDSDS+GCDIDS -GCGDODS 

+GCGSOGS -GCS IGS -GCDIGS -GCB IGS +GCGDOGS 

+GBD+GDSDS+GDSGS+GCGSOGS -GCS IDS -GCS IGS 
-GCBIDS-GCBIGS 

+GSUB 

+GBD+GBFEB+GCTEB -GC JEEB+GCDEEB 

+GDR+GBFEB+GBFCB+GBREB+GBRCB-GCJEEB 

-GCJCCB+GCDEEB+GCDECB+GCDCEB+GCDCCB 

+GSUB+GBRCB-GCTCB-GCJCCB+GCDCCB 

+GDSGS+GCDIGS -GCGDOGS 

-GDSDS -GDSGS -GCDIDS -GCDIGS+GCGDODS 
+GCGDOGS 

-GDR 

-GCSIDS-GCDIDS-GCBIDS+GCGDODS 

-GCGSOGS+GCSIDS+GCSIGS+GCDIDS+GCDIGS 
+GCB I DS +GCB IGS -GCGDODS -GCGDOGS 

-GDSDS+GCSIDS+GCBIDS 

-GDSGS -GCGSOGS+GCS IGS+GCB IGS 

-GBD 

-GSUB 

-GBD 

-GBFEB -GBFCB -GCTEB -GCTCB+GC JEEB -GCDEEB 
-GCDECB 
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VALUE ( LOCY ) =VALUE ( LOCY ) +GBFCB+GCTCB+GCDECB 
LOCY=LVN+NODPLC (LOC+25 ) 

VALUE { LOCY ) =VALUE ( LOCY ) -GDR 
LOCY=LVN+NODPLC (LOC+26 ) 

VALUE (LOCY) =VALUE (LOCY) -GBFEB-GBREB+GCJEEB-GCDEEB-GCDCEB 
LOCY=LVN+NODPLC (LOC+27 ) 

VALUE (LOCY) =VALUE (LOCY) -GBFCB-GBRCB+GCJCCB-GCDECB-GCDCCB 
LOCY=LVN+NODPLC (LOC+28) 

VALUE ( LOCY ) =VALUE ( LOCY ) -GSUB 
LOCY=LVN+NODPLC ( LOC+2  9 ) 

VALUE ( LOCY ) =VALUE ( LOCY ) +GBREB-GCTEB+GCDCEB 
LOCY=LVN+NODPLC ( LOC+3  0 ) 

VALUE (LOCY) =VALUE (LOCY) -GBREB-GBRCB+GCTEB+GCTCB+GCJCCB-GCDCEB 
Sc  -GCDCCB 

1000  LOC=NODPLC(LOC) 

GO  TO  10 
END 


C.2  Subroutine  LDTMOS.f 


SUBROUTINE  LDTMOS (VDS, VGS, LOC, IDS,QSI,QDI,QBI,QGDOL) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

C 

C THIS  ROUTINE  EVALUATES  THE  DRAIN-SOURCE  CURRENT  (IDS), 

C THE  INTRINSIC  GRADED -CHANNEL  CHARGES  (QSI,  QDI  & QBI) , AND 

C THE  GATE-DRAIN  OVERLAP  CHARGE  (QGDOL)  IN  THE  LATERAL  IS40S  TRANSISTOR. 

C 

C MODELED  BY  YEONBAE  CHUNG  7/15/93 

C PROGRAMMED  BY  YEONBAE  CHUNG  8/24/94 

C UPGRADED  (VERSION-1.1)  BY  YEONBAE  CHUNG  9/27/94 

C 

C SPICE  VERSION  2G.6  SCCSID=STATUS  3/15/83 

COMMON  /STATUS/  OMEGA, TIME, DELTA,DELOLD(7) ,AG(7) ,VT,XNI,EGFET, 

& XMU, SFACTR,MODE,MODEDC, ICALC,  INITF, METHOD, IORD,MAXORD, NONCON, 

Sc  ITERNO, ITEMNO, NOSOLV, MODAC, IPIV, IVMFLG, IPOSTP, ISCRCH, lOFILE 

C SPICE  VERSION  2G.6  SCCSID=KNSTNT  3/15/83 

COMMON  /KNSTNT/  TWOPI,XLOG2,XLOG10,ROOT2, RAD, BOLTZ, CHARGE, CTOK, 

Sc  Ca4IN,  RELTOL,  ABSTOL,  VNTOL,  TRTOL,  CHGTOL,  EPSO , EPSSIL,  EPSOX, 

Sc  PIVTOL,PIVREL 

C SPICE  VERSION  2G.6  SCCSID=BLANK  3/15/83 
COMMON  /BLANK/  VALUE (200000) 

C 

DOUBLE  PRECISION  NI, LCH, LDl, LD2 , LD3 , NAO, NDD, KFO , KFl , KF2 , KF3 , KF4 , 

Sc  N4 , LAMDAD,  NO , NABBAR,  KBI , ICHl , ICHLIN,  IDS , IDR,  ICH2 , ICHSAT, 

Sc  K1,L1,K2,L2,K3,L3,K4,L4 
INTEGER  I, J,MAXIT5, IFLAG 
INTEGER  NODPLC(64) 

COMPLEX  CVALUE(32) 

EQUIVALENCE  (VALUE ( 1 ) , NODPLC ( 1 ) , CVALUE ( 1 ) ) 

C 
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DIMENSION  Z41 (0:99) , 242(0:100),  243(0:100), 

& 25(1:201),  G5(l:200) 

G41(W1,W2,W3,W4,W5,W6,W7)  = Wl/(  (W3+W6*DTAN(W4) ) *W7  ) 

& - W2/ (W3+W6*DTAN(W4) ) + W5 

G42(W8)  = -W8 
C 
C 

LOCV=NODPLC ( LOC+ 1 ) 

LOCM=NODPLC (LOC+9 ) 

TYPE=NODPLC (LOCM+2 ) 

LOCM=NODPLC ( LOCM+ 1 ) 


Q-k  -k  -k  it  it  it  it  it  it  it  if  it  it  it  ir  ir  ir  it  ir  ir  it  ir  it  ir  ir  ir  it  it  ir  it  it  it  if  it  it  it  ir  ir  it  it  it  it  it  it  it  it  it  it  it  "k  it  it  it  it  it  -k  ir  it  1r  it  ir  ir  "k  ir  ir  it  it  ir  ir  irQ 

C PHYSICAL  CONSTANTS  C 

Qit  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  k k it Q 


Q 

= 

CHARGE 

NI 

= 

XNI 

EOX 

EPSOX 

ESI 

= 

EPSSIL 

VTO 

= 

VT 

PI 

— 

TWOPI/2.0D0 

GAMMA 

= 

PI/6.0d0 

^k  k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k kQ 

C MODEL  PARAMETERS  C 

Qk  k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k kQ 

wz 

= 

VALUE (LOCV+l) 

LCH 

= 

VALUE (LOCM+1) 

LDl 

= 

VALUE (LOCM+2) 

LD2 

= 

VALUE (LOCM+3 ) 

LD3 

= 

VALUE(LOCM+4) 

TOX 

= 

VALUE(LOCM+6) 

TFOX 

VALUE(LOCM+7) 

YJB 

= 

VALUE(LOCM+8) 

NAO 

VALUE (LOCM+ 11) 

NDD 

VALUE (LOCM+ 12) 

VSAT 

= 

VALUE (LOCM+ 15) 

THETA 

rr 

VALUE (LOCM+16) 

UNCHO 

= 

VALUE (LOCM+ 17) 

UND 

= 

VALUE (LOCM+ 18) 

BETAN 

VALUE (LOCM+ 19) 

VFBCH 

— 

VALUE (LOCM+2 6) 

VFBDR 

VALUE (LOCM+27) 

PSIBD 

= 

VALUE (LOCM+2 8) 

FQD 

VALUE (LOCM+2 9) 

KFO 

= 

VALUE (LOCM+3 0) 

KFl 

= 

VALUE (LOCM+3 1) 

KF2 

= 

VALUE (LOCM+3 2) 

KF4 

VALUE (LOCM+3 3) 

^k  k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k^ 

C ITERATION  CONSTANTS  C 

Qk  k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k kQ 

N4  = 30.0D0 
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MAXIT5  = 30 

EPS5  = l.OD-3 

Q ★★★★«>**★*★***#★**•*****<#*  i«r  *************  lAr  ****★★  lAr  lAr  lAr 

C CALCULATION  OF  PHYSICAL  QUANTITIES  C 

Q*************1t**1iii-kit*-IHr1t1c1i1r***1t*1tiHc1t*itit*1t1t***1t********************ir**-kQ 

COX  = EOX/TOX 

CFOX  = EOX/TFOX 

PIBO  = VTO*DLOG(NAO/NI) 

VTH  = VFBCH  + 2.0D0*PIB0  + DSQRT (4 . 0D0*Q*ESI*NA0*PIB0) /COX 
ETACH  = DLOG(NAO/NDD) 

UNCH  = UNCH0/(  l.ODO  + THETA* (VGS- VTH)  ) 

ECC  = VS AT/ UNCH 
IF  ( VGS  .LT.  1.1D0*VTH  ) THEN 
KF3  = 10.0D0*KF0*(VGS-VTH) /VTH 
ELSE 

KF3  = KFO 
ENDIF 

CDO  = KF3*DSQRT(  Q*ESI*NA0/ (4 . 0D0*PIB0)  ) 

VCHSAT  = KF4*(  COX* (VGS-VTH)  + CD0*ETACH*PIB0  )/(  COX  + CDO  ) 

ECD  = VSAT/UND 

IF  ( VDS  .GT.  -PSIBD  ) THEN 

YB  = DSQRT(  2.0D0*ESI*(PSIBD+VDS)/(Q*NDD)  ) 

ELSE 

YB  = O.ODO 
ENDIF 

A1  = (TFOX-TOX)/(2.0DO*LD2) 

B1  = YJB  + LDl 

D1  = (TFOX-TOX)/(  2.0D0*(LD1+LD2)  ) 

RC2  = DSQRT(  (YJB+LD1+LD2) **2  + 0 . 25D0* (TFOX-TOX) **2  ) 

XI  = YJB  + LDl  + LD2 

X3  = YJB  + LDl  + LD2  + LD3 

LAMDAD  = DSQRT(  ESI*VT0/ (Q*NDD)  ) 

HO  = KF1*LAMDAD 

NO  = Q*NDD/ESI 

PIBLCH  = VTO*DLOG(NDD/NI) 

PIBBAR  = 0.5D0*(  PIBO  + PIBLCH  ) 

NABBAR  = DSQRT ( NA0*NDD  ) 

KBI  = (2.0D0/3.0D0)*DSQRT(2,0D0*Q*ESI*NABBAR) 

ALPHA  = DATAN2(  (TFOX-TOX) , (2 . 0D0*LD2)  ) 


it  if  it  it  it  it  it  it  it  it  1r  ft  it  ir  it  ir  it  it  ir  ir  it  -k  it  ir  "k  it  It  it  it  it  ir  it  It  it  1r  it  it  it  it  It  it  it  it  ft  it  it  it  it  it  ir  it  ir  ir  ir  ir  It  it  "k  it  "k  ir  if  ir  ir  it  ir  "k  ir  itQ 

C COMPUTE  IDS  C 

Q ********************#***********************************'***********'** 


C..FOR  THE  CASE  OF  VDS  <=  0.0  OR  VGS  <=  VTH 

IF  ( ( VDS  .LE.  O.ODO  ) .OR.  ( VGS  .LE.  VTH  ) ) THEN 
IDS  = O.ODO 
GO  TO  300 
ENDIF 

C ==  EVALUATE  IDS  WHEN  VDS  > 0.0  & VGS  > VTH  BY  SECANT  METHOD  ====  C 


Z5(l)  = 0.1D0*VDS 
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Z5(2)  = 0.9D0*VDS 
DO  230  J = 1,  MAXIT5 
VCH  = VDS  - Z5(J) 

C . . CALCULATION  OF  CHANNEL  CURRENT 
IFLAG  = 0 

IF  ( VCH  .GT.  O.ODO  ) THEN 
IF  ( VCH  .LT.  VCHSAT  ) THEN 

UNLIN  = UNCH/(  1.0D0+ (VCH/ (ECC*LCH) ) **BETAN  ) ** (l.ODO/BETAN) 
ICHl  = 2.0D0*COX*(VGS-VTH)  + 2 . 0D0*ETACH*CD0*PIB0 
& - (COX+CDO) *VCH 

ICHLIN  = 0.5D0*WZ*UNLIN*ICH1*VCH/LCH 
IDS  = ICHLIN 
IDR  = ICHLIN 
ELSE 

UNSAT  = UNCH 

& / (1 . 0D0+ (VCHSAT/ (ECC*LCH) ) **BETAN) ** ( 1 . ODO/BETAN) 

ICH2  = ( COX*(VGS-VTH)  + CD0*ETACH*PIB0  )**2  / ( COX+CDO  ) 
ICHSAT  = 0.5D0*WZ*UNSAT*KF4*(2.0D0-KF4)*ICH2/LCH 
IDS  = ICHSAT 
IDR  = ICHSAT 
ENDIF 
ELSE 

UNLIN  = UNCH* (l.ODO -VCH) 

ICHl  = 2.0D0*COX*(VGS-VTH)  + 2 . 0D0*ETACH*CD0*PIB0 
& - (COX+CDO)* VCH 

ICHLIN  = 0.5D0*WZ*UNLIN*ICH1*VCH/LCH 
IDS  = ICHLIN 
IDR  = ICHLIN 
IFLAG  = 1 
ENDIF 

C . . CALCULATION  OF  EU ( 0 ) 

IF  ( VCH  .GT.  O.ODO  ) THEN 
VUO  = UNCH* (KF2* VCH /LCH) 

& /(  l.ODO  + ( KF2*VCH/(ECC*LCH)  )**BETAN  )**( 1 . ODO/BETAN) 

ELSE 

VUO  = UNCH* ( l.ODO -VCH )*KF2*VCH/LCH 
ENDIF 

EUO  = - VU0*ECD/(  VSAT  - VUO  ) 

C.. FINDING  Q-POINT 

IF  ( VDS  .LT.  (VGS-VFBDR)  ) THEN 
IF  ( YB  .LE.  LDl  ) THEN 
XQ  = YJB  + YB 
YQ  = O.ODO 

ELSEIF  ( ( LDl  .LT.  YB  ) .AND.  ( YB  .LE.  (RC2-YJB)  ) ) THEN 
XQl  = DSQRT(  (1.0D0+A1**2)*(YJB+YB)**2  - (A1*B1)**2  ) 

XQ  = ( A1*A1*B1  + XQl  ) / ( l.ODO  + A1*A1  ) 

YQ  = A1*(XQ-B1) 

ELSEIF  ( YB  .GT.  ( RC2-YJB  ) ) THEN 

XQ2  = DSQRT(  ( 1 . 0D0+D1**2 ) * (YJB+YB) **2  - (D1*YJB)**2  ) 

XQ  = ( D1*D1*YJB  + XQ2  ) / ( l.ODO  + D1*D1  ) 

YQ  = D1*(XQ-YJB) 
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Sc 

Sc 

Sc 


Sc 

110 

Sic 

Sic 

120 

Sic 


130 


ENDIF 

ELSE 

YOX  = DSQRT(  (ESI/COX) **2 

+ 2.0D0*ESI*(VDS+VFBDR-VGS)/(Q*NDD)  ) - ESI/COX 
YFOXl  = DSQRT(  (ESI/CFOX) **2 

+ 2.0D0*ESI*(VDS+VFBDR-VGS)/(Q*NDD)  ) 

YFOX2  = DSQRT(  (ESI/CFOX) **2  - (ESI/COX) **2 

+ ( ESI/COX  + 0.5DO*(TFOX-TOX)  )**2  ) 

YFOX  = YFOXl  - YFOX2 

Cl  = ( YOX  - YFOX  - 0.5DO*(TFOX-TOX)  ) / LD2 

X2  = YJB  + LDl  + LD2  + 2 . ODO* (LD1+LD2 ) *YFOX/ (TFOX-TOX) 

IF  ( YOX  .LE.  0.5DO*(TFOX-TOX)  ) THEN 
IF  ( ( YJB  + YB  ) .GE.  YOX  ) THEN 

XQ3  = DSQRT(  (YJB+YB)**2  - YOX*YOX  ) 

ELSE 

GO  TO  110 
ENDIF 

IF  ( ( (YJB+YB)  .LE.  RC2  ) .AND.  ( YJB  .LT.  XQ3  ) 

.AND.  ( XQ3  .LE.  (YOX/Al+Bl)  ) ) THEN 
XQ  = XQ3 
YQ  = YOX 
GO  TO  210 
ENDIF 

IF  ( ( DSQRT(1.0D0+A1*A1) * (YJB+YB)  ) .GE.  (A1*B1)  ) THEN 
XQ4  = ( A1*A1*B1  + DSQRT(  ( 1 . 0D0+A1*A1) * (YJB+YB) **2 

- A1*A1*B1*B1  ) ) / ( l.ODO  + A1*A1  ) 

ELSE 

GO  TO  120 
ENDIF 

IF  ( ( (YJB+YB)  .LE.  RC2  ) .AND. 

( XQ4  .GT.  (YOX/Al+Bl)  ) ) THEN 
XQ  = XQ4 
YQ  = A1*(XQ4-B1) 

GO  TO  210 
ENDIF 

IF  ( ( DSQRT(1.0D0+D1*D1)* (YJB+YB)  ) .GE.  (D1*YJB)  ) THEN 
XQ5  = ( D1*D1*YJB  + DSQRT(  (1.0D0+D1*D1) * (YJB+YB) **2 

- (D1*YJB)**2  ) ) / ( l.ODO  + D1*D1  ) 

ELSE 

GO  TO  130 
ENDIF 

IF  ( (YJB+YB)  .GT.  RC2  ) THEN 
XQ  = XQ5 

YQ  = D1*(XQ5-YJB) 

GO  TO  210 
ENDIF 
GO  TO  240 
ELSE 

IF  ( (YJB+YB)  .GE.  YOX  ) THEN 

XQ6  = DSQRT(  (YJB+YB) **2  - YOX*YOX  ) 

ELSE 

GO  TO  140 
ENDIF 

IF  ( ( YJB  .LT.  XQ6  ) .AND.  ( XQ6  .LE.  (YJB+LDl)  ) ) THEN 
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140 

Sc 

Sc 

150 


160 

Sc 


170 


180 


XQ  = XQ6 
YQ  = YOX 
GO  TO  210 
ENDIF 

IF  ( (DSQRT(1.0D0+C1*C1)*(YJB+YB))  .GE.  (Y0X+B1*C1)  ) THEN 
XQ7  = ( C1*(Y0X+B1*C1)  + DSQRT( (1.0D0+C1*C1)*(YJB+YB)**2 
- (Y0X+B1*C1) **2)  ) / ( l.ODO  + C1*C1  ) 

ELSE 

GO  TO  150 
ENDIF 

IF  ( ( (YJB+LDl)  .LT,  XQ7  ) .AND. 

( XQ7  .LE.  XI  ) ) THEN 
XQ  = XQ7 

YQ  = YOX  - C1*(XQ7-B1) 

GO  TO  210 
ENDIF 

IF  ( (XI  .LT.  X2  ) .AND.  ( X2  .LE.  X3  ) ) THEN 

IF  ( (YJB+YB)  .GE.  ( YFOX  + 0 . 5D0* (TFOX-TOX)  ) ) THEN 
XQ8  = DSQRT( (YJB+YB) **2-{YFOX+0.5D0* (TFOX-TOX) ) **2) 
ELSE 

GO  TO  160 
ENDIF 

IF  ( (XI  .LT.  XQ8  ) .AND.  ( XQ8  .LE.  X2  ) ) THEN 
XQ  = XQ8 

YQ  = YFOX  + 0.5D0* (TFOX-TOX) 

GO  TO  210 
ENDIF 

IF  ( (DSQRT(1.0D0+D1*D1) * (YJB+YB) ) .GE.  (D1*YJB)  ) THEN 
XQ9  = ( D1*D1*YJB  + DSQRT(  (1 . 0D0+D1*D1) * (YJB+YB) **2 
- (D1*YJB)**2  ) ) / ( l.ODO  + D1*D1  ) 

ELSE 

GO  TO  170 
ENDIF 

IF  ( XQ9  .GT.  X2  ) THEN 
XQ  = XQ9 

YQ  = D1*(XQ9-YJB) 

GO  TO  210 
ENDIF 
ENDIF 

IF  ( X2  .GT.  X3  ) THEN 

IF  ( (YJB+YB)  .GE.  ( YFOX  + 0 . 5D0* (TFOX-TOX)  ) ) THEN 
XQIO  = DSQRT((YJB+YB)**2-(YFOX+0.5DO*(TFOX-TOX))**2) 
ELSE 

GO  TO  180 
ENDIF 

IF  ( (XI  .LT.  XQIO  ) .AND.  ( XQIO  .LE.  X3  ) ) THEN 
XQ  = XQIO 

YQ  = YFOX  + 0.5D0* (TFOX-TOX) 

GO  TO  210 
ENDIF 
ENDIF 
GO  TO  240 
ENDIF 
ENDIF 
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C . . CALCULATION  OF  VDR  BY  RUNGE-KUTTA  METHOD 


210 


220 


WSCR  = DSQRT(  (XQ-YJB) **2  + YQ*YQ  ) 
AO  = IDR/ (2.0D0*ESI*UND*WZ) 

BO  = IDR/(2.0D0*ESI*ECD*UND*WZ) 

H4  = WSCR/N4 
Z42(0)  = EUO 
Z43 ( 0 ) = VCH 

DO  220  1=0,  IDINT(N4)  - 1 
Z41(I)  = FLOAT(I)*H4 


K1  = H4*G41( 

AO, 

BO, 

HO, 

GAMMA, 

NO 

LI  = H4*G42( 

Z42(I) 

) 

K2  = H4*G41( 

AO, 

BO, 

HO, 

GAMMA, 

NO 

Z41(I),  Z42(I)  ) 


Z4KD+H4/2.0D0,  Z42(I)+K1/2.0D0  ) 

L2  = H4*G42(  Z42 ( I) +K1/2 . ODO  ) 

K3  = H4*G41(  AO,  BO,  HO,  GAMMA,  NO, 

Z41(I)+H4/2.0D0,  Z42(I)+K2/2.0D0  ) 

L3  = H4*G42(  Z42 ( I) +K2/2 . ODO  ) 

K4  = H4*G41(  AO,  BO,  HO,  GAMMA,  NO,  Z41(I)+H4,  Z42(I)+K3  ) 
L4  = H4*G42(  Z42(I)+K3  ) 

Z42{I+1)  = Z42(I)  + (Kl+2.ODO*K2+2.0DO*K3+K4)/6.0DO 
Z43(I+1)  = Z43(I)  + (L1+2.0D0*L2+2.0D0*L3+L4) /6.0D0 
CONTINUE 

VSCR  = Z43(N4)  - VCH 
VDRP  = VSCR 


C 

C 

G5(J)  = Z5(J)  - VDRP 

IF  ( DABS(  Z5(J)  - VDRP  ) .LT.  EPS5  ) THEN 
VDR  = Z5(J) 

GO  TO  240 
ENDIF 

IF  ( J .EQ.  1 ) THEN 
GO  TO  230 
ENDIF 

Z5(J+1)  = Z5(J)  - ( Z5(J)  - Z5(J-1)  ) * G5(J) 
Sc  / ( G5(J)  - G5(J-1)  ) 

230  CONTINUE 
GO  TO  240 


C ==  END  OF  DETERMINING  IDS  FOR  NORMAL  BIAS  CONDITIONS  ===========  C 

C . . CHECK  THE  SOLUTION 

240  IF  ( IFLAG  .EQ,  1 ) THEN 
IDS  = O.ODO 
ENDIF 

p ********************************«**********************«************'* 

C COMPUTE  QSI  C 

p**********************************************************************Q 

300  IF  { VGS  .LE.  VTH  ) THEN 
QSI  = O.ODO 
ELSE 

IF  ( (VDS.LE. O.ODO)  .OR.  (VCH .LE. 0 . ODO)  ) THEN 
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QSI  = -(WZ*LCH/6.0D0) 

& *(  3.0D0*COX*(VGS-VTH)  + 2 . ODO*CDO*ETACH*PIBO  ) 

ELSEIF  ( VCH  .LT,  VCHSAT  ) THEN 

QSI  = -(WZ*LCH/6.0D0)*(  3 . ODO*COX* (VGS-VTH) 

& + 2.0D0*CD0*ETACH*PIB0  - (COX+CDO) *VCH  ) 

ELSEIF  ( VCH  .GE.  VCHSAT  ) THEN 
QSI  = -(WZ*LCH/6.0D0) 

& * ( ( 3 . ODO -KF4 ) *COX* ( VGS-VTH ) + ( 2 . ODO -KF4 ) *CD0  *ETACH*PIBO ) 

ENDIF 
ENDIF 

C COMPUTE  QDI  C 

IF  ( VGS  .LE.  VTH  ) THEN 
QDI=  O.ODO 
ELSE 

IF  ( (VDS.LE. O.ODO)  .OR.  (VCH .LE. 0 . ODO ) ) THEN 
QDI  = -(WZ*LCH/6.0D0) 

St  *(  3 .ODO *COX*{ VGS-VTH)  + 4 . 0D0*CD0*ETACH*PIB0  ) 

ELSEIF  { VCH  .LT.  VCHSAT  ) THEN 

QDI  = -(WZ*LCH/6.0D0)*(  3 . ODO *COX* (VGS-VTH) 

& + 4.0D0*CD0*ETACH*PIB0  - 2 . ODO* (COX+CDO) *VCH  ) 

ELSEIF  ( VCH  .GE.  VCHSAT  ) THEN 
QDI  = -(WZ*LCH/6.0D0) 

St  *(  (3 .0D0-2.0D0*KF4)  *COX*  (VGS-VTH) 

St  + (4.0D0-2.0D0*KF4)  *CD0*ETACH*PIB0  ) 

ENDIF 

ENDIF 


Qir  ir  ir  ir  ir  it  it  it  it  ir  ir  it  it  ir  ir  it  it  it  it  it  it  it  it  it  it  ir  ir  ir  it  it  it  it  it  it  it  it  it  ir  it  it  it  ir  it  it  it  it  -it  ir  ir  ir  "k  ir  it  it  ir  it  ir  it  it  it  ir  it  it  ir  it  it  it  it  ir  icQ 

C COMPUTE  QBI  C 

Qit  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it Q 

IF  ( VGS  .LE.  VFBCH  ) THEN 
QBI  = WZ*LCH*COX*(VFBCH-VGS) 

ELSEIF  ( ( VFBCH  .LT.  VGS  ) .AND.  ( VGS  .LE.  VTH  ) ) THEN 
QBI  = -1.5D0*WZ*LCH*KBI*DSQRT(  2.0D0*PIBBAR  ) 

& *DSQRT(  (VGS-VFBCH)/( VTH -VFBCH)  ) 

ELSEIF  ( VGS  .GT.  VTH  ) THEN 

IF  ( (VDS.LE. O.ODO)  .OR.  (VCH .LE. 0 . ODO)  ) THEN 
QBI  = -1.5D0*WZ*LCH*KBI*DSQRT(  2.0D0*PIBBAR  ) 

ELSEIF  ( VCH  .LT.  VCHSAT  ) THEN 
QBI  = -WZ*LCH*KBI 

& * ( (VCH+2 . 0D0*PIBBAR) **1 . 5- (2 . 0D0*PIBBAR) **1 . 5) /VCH 

ELSEIF  ( VCH  .GE.  VCHSAT  ) THEN 

QBIl  = ( (COX+CDO) /KF4)  / ( COX* (VGS-VTH)  + CD0*ETACH*PIB0  ) 
QBI2  = KF4*(  COX* (VGS-VTH)  + CD0*ETACH*PIB0  )/ (COX+CDO) 

& + 2.0D0*PIBBAR 

QBI  = -WZ*LCH*KBI*QBI1*(  QBI2**1.5  - (2 . 0D0*PIBBAR) **1 . 5 ) 
ENDIF 
ENDIF 


Qit  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it Q 

C COMPUTE  QGDOL  C 
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^**********************************************************************p 

IF  ( VDS  .LT,  (VGS-VFBDR)  ) THEN 

QGDOL  = -WZ*(  COX*LDl  + 0 . 5D0* (COX+CFOX) *LD2  + CFOX*LD3  ) 

& *l  VDS  - VGS  + VFBDR  ) 

ELSE 

YOX  = DSQRT(  (ESI/COX) **2 

& + 2.0D0*ESI*(VDS+VFBDR-VGS)/(Q*NDD)  ) - ESI/COX 

YFOXl  = DSQRT(  (ESI/CFOX) **2 

& + 2.0D0*ESI*(VDS+VFBDR-VGS)/(Q*NDD)  ) 

YFOX2  = DSQRT(  (ESI/CFOX) **2  - (ESI/COX) **2 

& + ( ESI/COX  + 0.5DO*(TFOX-TOX)  )**2  ) 

YFOX  = YFOXl  - YFOX2 

IF  ( YOX  .LE.  0.5DO*(TFOX-TOX)  ) THEN 

QGDOL  = -Q*FQD*WZ*NDD*YOX* ( LDl  + LD2*YOX/ (TFOX-TOX)  ) 

ELSE 

QGDOL  = -Q*FQD*WZ*NDD 

& *(  (LD1+0.5D0*LD2) *YOX  + (0 . 5D0*LD2+LD3 ) *YFOX  ) 

ENDIF 

ENDIF 

C 

C 

RETURN 

END 


C . 3 Subrout ine  LDTB JT . f 


SUBROUTINE  LDTBJT(VEB, VCB, LOC, IBF, IBR, ICT, QDE, QDC) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

C 

C THIS  ROUTINE  EVALUATES  THE  PARASITIC  BJT  CURRENTS  (IBF,  IBR 

C Sc  ICT),  AND  THE  FORWARD-BIASED  INJECTION  CHARGES  (QDE  Sc  QDC) 

C IN  THE  LATERAL  EMOS  TRANSISTOR. 

C 

C MODELED  BY  YEONBAE  CHUNG  7/15/93 

C PROGRAMMED  BY  YEONBAE  CHUNG  8/24/94 

C 

C SPICE  VERSION  2G.6  SCCSID=STATUS  3/15/83 

COMMON  / STATUS / OMEGA , TIME , DELTA , DELOLD (7),AG(7),VT,XNI, EGFET , 

Sc  XMU, SFACTR,MODE,MODEDC, ICALC, INITF, METHOD, IORD,MAXORD, NONCON, 

Sc  ITERNO, ITEMNO, NOSOLV, MODAC , IPIV, IVMFLG, IPOSTP, ISCRCH, lOFILE 
C SPICE  VERSION  2G.6  SCCSID=KNSTNT  3/15/83 

COMMON  /KNSTNT/  TWOPI, XLOG2 , XLOGIO , ROOT2 , RAD, BOLTZ, CHARGE, CTOK, 

& GMIN, RELTOL, ABSTOL, VNTOL, TRTOL, CHGTOL, EPSO , EPSSIL, EPSOX, 

& PIVTOL, PIVREL 

C SPICE  VERSION  2G.6  SCCSID=BLANK  3/15/83 
COMMON  /BLANK/  VALUE (200000) 

C 

DOUBLE  PRECISION  NI,NDD,NAS,NABEFF,LAB,LAC, JNOYSE, JNOYSC, JRCBO, 

& ICTl , ICT2 , ICT3 , ICT, IBFl , IBF2 , IBF,  IBRl , IBR2 , IBR 
INTEGER  NODPLC(64) 
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COMPLEX  CVALUE(32) 

EQUIVALENCE  (VALUE ( 1 ) , NODPLC ( 1 ) , CVALUE ( 1 ) ) 

C 

C 

LOCV=NODPLC ( LOC+ 1 ) 

LOCM=NODPLC (LOC+9 ) 

TYPE=NODPLC (LOCM+2 ) 

LOCM=NODPLC { LOCM+ 1 ) 

p**********************************************************************C 

C PHYSICAL  CONSTANTS  C 

p**********************************************************************c 

Q = CHARGE 

NI  = XNI 

ESI  = EPSSIL 

VTO  = VT 


C MODEL  PARAMETERS  C 


AE 

VALUE(LOCV+8) 

AC 

= 

VALUE (LOCV+9 ) 

AEC 

VALUE (LOCV+ 10) 

YJDP 

VALUE(LOCM+9) 

YEPI 

= 

VALUE (LOCM+10) 

NDD 

VALUE (LOCM+ 12 ) 

NAS 

= 

VALUE (LOCM+ 13) 

NABEFF 

= 

VALUE (LOCM+ 14) 

BB 

VALUE ( LOCM+2 0) 

TAUNE 

= 

VALUE (LOCM+2 1) 

TAUHB 

VALUE (LOCM+2 2 ) 

TAUHC 

= 

VALUE (LOCM+2 3) 

LAB 

= 

VALUE (LOCM+2 4) 

LAC 

= 

VALUE (LOCM+2 5) 

PSIBD 

VALUE (LOCM+28) 

Q* ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ 

C CALCULATION  OF  PHYSICAL  QUANTITIES  C 

p**********************************************************************C 


PSISUB  = VTO*DLOG(NDD*NAS/ (NI*NI) ) 
JNOYSE  = Q*NI*NI*YJDP/ (NABEFF*TAUNE) 
JNOYSC  = -Q*NI*LAC/TAUHC 
JRCBO  = Q*NI*LAB/TAUHB 


IF  ( VEB  .LT. 

YE  = YJDP  + 
ELSE 

YE  = YJDP 
ENDIF 

IF  ( VCB  .LT. 
YC  = YEPI  - 


PSIBD 
DSQRT ( 


) THEN 

2 . 0D0*ESI* (PSIBD-VEB) / (Q*NDD) 


PSISUB 
DSQRT ( 


) THEN 

2 . 0D0*ESI*NAS* (PSISUB-VCB) 
/ (Q*NDD* (NAS+NDD) ) ) 


ELSE 
YC  = 
ENDIF 


) 


YEPI 
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WBN  = (YC-YE)/LAB 


(;************************«*****************************************^**^^^ 

C PREVENT  UNDERFLOW  C 

C**********************************************************************p 

IF  ( VEB/VTO  .LT.  -200. ODO  ) THEN 
EPVEBl  = O.ODO 
ELSE 

EPVEBl  = DEXP (VEB/VTO ) 

ENDIF 

IF  ( VEB/(2.0D0*VT0)  .LT.  -100. ODO  ) THEN 
EPVEB2  = O.ODO 
ELSE 

EPVEB2  = DEXP ( VEB/ { 2 . ODO  *VT0 ) ) 

ENDIF 

IF  ( VCB/ ( 2 . ODO *VT0 ) .LT.  -100. ODO  ) THEN 
EPVCB2  = O.ODO 
ELSE 

EPVCB2  = DEXP(VCB/(2.0D0*VT0)) 

ENDIF 


Q ********  «r  **#*****«******************«***#*******************'********« 

C COMPUTE  IBF  C 

Q**-k**ft’kirir'kiritir-kiritiriritiriririritititititititir’kirir-kititititiritiririr’k'kit’kirir1riritiririririririr**itit*irir’kiritQ 

IF  ( WBN  .GT.  O.ODO  ) THEN 

IBFl  = AE*JRCB0*(  DCOSH(WBN)  - 1 . ODO  )/{  DSINH(WBN)  ) 

Sc  * ( EPVEB2  - l.ODO  ) 

ELSE 

IBFl  = O.ODO 
ENDIF 

IBF2  = AE*JN0YSE*(  EPVEBl  - l.ODO  ) 

IBF  = IBFl  + IBF2 


^**********************************************************************P 

C COMPUTE  IBR  C 

Qir**it*ir’kirftir*iririririrititirititititirititiritirit**itititic*ititirit*itir1r’k-kititiritiririrititirititiriritirititiririticirQ 

IF  { WBN  .GT.  O.ODO  ) THEN 

IBRl  = JRCBO* (DCOSH (WBN) -l.ODO) /DSINH (WBN)  - JNOYSC 
ELSE 

IBRl  = -JNOYSC 
ENDIF 

IBR2  = EPVCB2  - l.ODO 
IBR  = AC*IBR1*IBR2 

Qir  ir  ir  it  ir  it  ir  ir  it  ir  it  it  ir  it  it  ir  it  it  ir  It  it  it  it  it  -k  'k  1r  * it  it  ir  -k  -k  ir  it  it  ir  it  ir  -k 

C COMPUTE  ICT 

Q*  if  k it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  k it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it 

ICTl  = (AEC/BB) *JN0YSE*(  EPVEBl  - 1 
ICT2  = (AEC/BB) *JN0YSC* ( EPVCB2  - 1 
IF  ( WBN  .GT.  l.OD-3  ) THEN 

ICT3  = AEC*JRCB0*(  l.ODO  + DCOSH (WBN) /BB  )/(  DSINH(WBN)  ) 

St  *(  EPVEB2  - EPVCB2  ) 

ELSE 

ICT3  = AEC*JRCB0*(  1 . ODO+DCOSH ( 1 . OD-3 ) /BB  ) / (DSINH ( 1 . OD-3 ) ) 


kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkQ 

c 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkk^ 

.ODO  ) 

.ODO  ) 
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Jc  *{  EPVEB2  - EPVCB2  ) 

ENDIF 

ICT  = ICTl  + ICT2  + ICT3 

C COMPUTE  QDE  C 

Q-k*’k'k’k’k’k’k'k’ir'k'k'kir’kir'k’kit1r*ititiritit**irit**iritir’kit'k’kiritirir'k'kir'kiririritifkir’kiriHrit’kit’k*'k*ifk*ifkQ 

IF  ( WBN  .GT.  O.ODO  ) THEN 

QDEl  = AE*TAUHB*JRCBO* ( DCOSH(WBN)  - l.ODO  )/(  DSINH(WBN)  ) 

& * ( EPVEB2  - l.ODO  ) 

ELSE 

QDEl  = O.ODO 
ENDIF 

QDE2  = AE*TAUNE* JNOYSE* { EPVEBl  - l.ODO  ) 

QDE  = QDEl  + QDE2 


Q ************  «t  ******************************************************** 

C COMPUTE  QDC  C 

it  1r  it  it  it  it  It  it  It  -k  It  it  it  it  It  it  ir  ir  ir  it  it  it  it  it  it  ir  it  ir  it  it  it  it  it  ir  it  ir  ir  it  it  ir  it  it  it  ir  ir  it  it  it  it  it  ir  it  it  it  it  it  ir  -k  "k  it  ir  it  it  ir  ir 

IF  ( WBN  .GT.  O.ODO  ) THEN 

QDCl  = TAUHB*JRCB0*(DCOSH(WBN)-1.0D0)/DSINH(WBN) 

& - TAUHC*JNOYSC 

ELSE 

QDCl  = -TAtJHC*JNOYSC 
ENDIF 

QDC2  = EPVCB2  - l.ODO 
QDC  = AC*QDC1*QDC2 
C 
C 

RETURN 

END 
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